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1.0  SCOPE 

1.1  Identification 

The  Computer  Software  Configuration  Item  (CSCI),  identified  as  the  Globally  Relocatable  Navy 
Tide/Atmosphere  Modeling  System  (PCTides)  consists  of  a  2-  and  3-dimensional  barotropic 
tide/surge  model,  called  the  Global  Environmental  Modeling  Services  (GEMS)  Coastal  Ocean 
Model  (GCOM2D  and  GCOM3D),  and  a  Mesoscale  Atmospheric  Prediction  System  (MAPS). 

GCOM2D  is  a  depth-integrated  shallow  water  model  designed  to  characterize  sea  level  and  currents 
on  or  near  continental  shelves.  It  features  a  wetting  and  draining  algorithm  for  simulating  coastal 
flooding  due  to  tides  or  storm  surge. 

GCOM3D  is  the  three-dimensional  counterpart  to  GCOM2D.  It  is  a  barotropic  model  for 
applications  where  current  structure  with  vertical  depth  is  required  and  tidal  and  wind  forcing  are 
dominant.  Atmospheric  forcing  for  GCOM2D/3D  is  provided  by  an  existing  operational  Navy 

model,  by  the  MAPS  system,  by  an  analytical  hurricane  vortex  model,  or  by  direct  point 
observation. 

MAPS  is  a  hydrostatic  primitive  equations  model  designed  to  provide  high-resolution 
representations  of  anemometer  level  winds  and  surface  pressure  as  atmospheric  boundary 
conditions  for  GCOM2D  and  GCOM3D.  To  this  end,  the  turbulence  closure  scheme  has  been 
desired  to  allow  the  model  to  be  run  with  its  lowest  model  level  at  anemometer  height,  thus 
providing  a  direct  simulation  of  the  winds  at  this  level.  The  equations  of  motion  are  coded  in 
advective  form  and  solved  using  a  semi-implicit  time  differencing  scheme  ensuring  that  the  model 
is  both  robust  and  economical  to  run,  even  in  regions  of  steep  terrain. 

1.2  Document  Overview 

The  purpose  of  this  Software  Design  Description  (SDD)  is  to  describe  the  software  design  and 
code  of  the  Globally  Relocatable  Navy  Tide/Atmosphere  Modeling  System  (PCTides).  It 
includes  the  mathematical  formulation  and  solution  procedures  for  GCOM2D/3D  and  MAPS  as 
well  as  flow  charts  and  descriptions  of  the  programs,  subprograms,  and  common  blocks. 

PCTides  is  part  of  the  Oceanographic  and  Atmospheric  Master  Library  (OAML),  and  is  actively 
configuration  managed  under  the  direction  of  that  authority.  This  document,  along  with  the 
Software  Requirements  Specification  and  the  Software  Test  Description  form  the  standard 
documentation  package  for  the  OAML  PCTides. 
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3.0  CSCI-WIDE  DESIGN  DECISION 

The  MAPS  system  uses  a  direct  access  file  format  where  all  model  variables  at  all  levels  are 
stored  at  each  given  time.  Refer  to  section  5.4a. 

4.0  CSCI  ARCHITECTURAL  DESIGN 

4.1  CSCI  Components 

a)  PCTides  can  be  divided  into  three  main  software  units:  CSC  MAPS,  CSC  GCOM2D  and 
CSC  GCOM3D.  These  are  briefly  described  below  along  with  commonly  used  library 
subroutines  and  data  libraries  required  for  smooth  operation  of  PCTides. 

b,  c)  CSC  MAPS-  This  program  consists  of  23  subroutines  and  6  common  blocks,  the  bulk  of 
which  are  used  to  provide  high  resolution  spatial  representation  of  the  wind  Lid  pressure 
field  for  the  GCOM  ocean  models. 

CSC  GCOM2D-  This  program  consists  of  21  subroutines  and  6  common  blocks  that  aid  in 
modeling  sea  surface  height  and  mean  barotropic  current  structure  in  coastal  regions  due  to 
tides  and  atmospheric  forcing. 

CSC  GCOM3D  -  This  three-dimensional  z-coordinate  hydrod>'namic  model  consists  of  two 
subroutines  and  two  common  blocks  that  are  used  for  coastal  ocean  modeling  where 
barotropic  current  structure  with  depth  is  required.  The  remaining  components  of  the  model 
^e  GCOM2D  subroutines  and  common  blocks.  This  CSC  can  be  run  in  either  barotropic 
(no  thermal  and  density  variation)  or  baroclinic  modes  (temperature  and  salinity  are  model 
variables),  but  only  the  barotropic  mode  is  supplied  with  this  installation. 

Standard  Library  Routines  for  PCTides  programs: 

Input/output  routines  -The  seventeen  subroutines  fall  into  four  file  formats: 

a)  Direct  access  files-all  model  variables  at  all  levels  are  stored  at  each  given  time 
(9  subroutines). 

b)  Sequential  files-  used  for  storage  of  output  data  (4  subroutines). 

c)  Time  scries  files-  used  for  storing  station  infonnation  in  both  tlie  atmospheric 
and  ocean  models  (2  subroutines). 

d)  Topography  files-standard  ASCII  format  used  to  write  tltc  topography  and 
bathymetry  of  a  model  grid  (2  subroutines). 

Proiection  routines-  In  the  PCTides  system,  both  MAPS  and  GCOM2D/3D  utilize  a 
standard  Cartesian  grid  projection  system  for  tlic  solution  of  tlie  equations  of  motion.  Four 
subroutines  aid  in  producing  this  grid  system. 
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d) 


Date^ime  Routines-  Four  subroutines  are  used  to  convert  between  date  and  time,  which  are 
specified  on  input  files,  and  an  invariant  reference  time  used  by  the  model. 


Filtering  routines-  Two  subroutines  are  included  in  the  model  to  reduce  numerical  errors  and 
computational  instabilities  that  may  lead  to  a  spurious  growth  of  short  waves  that  may  mask 
the  actual  numerical  simulation  or,  in  severe  cases,  could  cause  catastrophic  instability. 


Miscellaneous  MAPS  routines-  Four  subroutines 
pertain  to  the  solution  of  semi-implicit  equations. 


in  the  atmospheric  prediction  model 


The  PCTides  model  is  relatively  new  and  has  not  yet  been  incorporated  into  the  OAML 
software  collection. 


e)  PCTides  takes  up  about  360  MB  of  disk  space  in  the  PC  Windows  environment  The 
program  requires  the  same  amount  of  disk  space  for  UNIX  operating  systems. 

4.2  Concept  of  Execution 


PCTides  IS  a  platform  independent  system  that  may  be  run  in  the  PC  Windows  environment 
executed  through  an  interactive  menu  or  at  the  command  prompt  under  DOS  mode  or  the  UNIX 

menu  allows  for  a  logical  structure  by  which  the  user  can 
establish  the  required  environment  for  carrying  out  a  simulation  with  either  GCOM2D  or  GCOM3D 
m  any  part  of  the  world’s  oceans.  Tides  and/or  winds  can  drive  PCTides, 

Execution  of  tWs  CSCI  begins  by  defining  the  model  domain  and  generating  a  grid  of  the  analysis 
region  A  bath>metiy  grid  is  then  created  and  tidal  boundary  conditions  are  derived  for  the  region 
frem  global  tidal  files  containing  the  amplitude  and  phase  of  the  tidal  constituent  in  time  zone  zero 
(Greenwich).  Wind  data  are  acquired  from  three  possible  sources;  1)  NOGAPS  or  some  hieher 
resolution,  2)  manual  data  entiy,  3)  hurricane  model  adaptations. 

Pnor  to  execution  of  the  GCOM  program,  several  parameters  must  be  specified  such  as  a  wind  flag. 

1  e  flag,  run  time  and  tidal  start  time,  and  output  file  time  step,  among  others.  In  addition,  station 

‘STATlONTnA^fm  successfiilly.  The 

S I  ATIONS.DAT  file  may  have  up  to  30  stations  defined  per  model  run. 

Next,  the  type  of  model  must  be  specified.  Both  GCOM2D  and  GCOM3D  may  be  run  using  the 

GCOM3D  if  near-surface  currents  were 
required.  GCOM2D  is  faster  and  is  used  to  predict  sea  levels  as  a  result  of  tidal  and  or  hurricane 

fCoding  simulate  coastal  region  inundation  due  to  tidal  ebb  and  storm  surge  or 

The  user  must  specify  output  display  options.  Display  options  allow  the  plotting  of  spatial  fields  or 
ime  senes  of  sea  levels  and  ocean  currents.  Output  stations,  at  which  time  series  prediction  of  sea 
levels  and  currents  arc  stored,  may  also  be  specified. 


6 


PSI  Technical  Report  SSC-005-00 


PCTides  SDD 


Assistance  for  the  system  may  be  obtained  in  the  PC  Windows  menu  where  the  menu  procedures 
are  also  supported  with  help  functions  for  all  actions  of  choice.  Help  is  accessed  by  pressing  the 
“FI”  key  while  the  cursor  is  on  the  input  box  in  question. 

A  diagram  illustrating  the  basic  logic  underlying  operation  of  the  CSCI  through  PC  Windows  or  the 
command  prompt  is  shown  in  (Figure  4.2-1). 


Main  Menu 


Sub  Menus 


I  GENERATE  GRID 
GET  BATHYMETRY 

GET  TIDES 

MAPS  WINDS 
NOGAPS  WINDS 

MANUAL  WINDS 

HURRICANE 

PARAMETERS  j 
STATIONS 

GCOM2D 
GCOM3D  1 


- ►  ASAGRID.EXE 

- ►  GCOMGRID.EXE 

- ►  TIDES.EXE 

- ►  DA2SEQ.EXE 

►  NRL2SEQ,EXE 
forecast! 

>  TRACK  -►  WINDS.EXE 


HINDCAST 

TRACK 


CYCLONE.EXE 


CYCLONE 

MODEL 


^^^-►CYCLONE.EXE 

GCOM.DAT 

STATIONS.DAT 

GCOM2D.EXE 

GCOM3D.EXE 


i 

! 

CURRENTS 

HEIGHTS 

SPEED  TS 

DIREC  TS 

TIDAL  BC 

2D  BATHYMETRY 

3D  BATHYMETRY 


OMDISP.EXE 

ZPLOT.EXE 

SPDPLOT.EXE 

DIRPLOT.EXE 

PLOTIDE.EXE 

BATHPLOT.EXE 

PLOT3D.EXE 


Figure  4.2-1:  Chart  illustratiim  the  PC  Windows  Menu  and  related  files. 
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4.3  INTERFACE  DESIGN 


4.3. 1  INTERFA CE IDENTIFICA  TION AND  DIA GRAMS 


The  only  Navy  standard  PCTides  external  interfaces  are  the  input  and  output  files. 
The  input  files  consist  of  the  following  types  of  data; 


1. 

2. 

3. 

4. 


archived  topographic  data, 
model  domain  parameters  in  ‘GRIDGEN.DAT’. 

Grenoble  tidal  data. 

wind  date  from  NOGAPS  or  some  available  higher  resolution,  from  manual  data  entry  or 
denved  from  the  hurricane  model. 

inputs  to  the  coastal  ocean  model  include  various  parameter  and  station  location  data. 


The  output  files  consist  of  the  following  types  of  data: 

1.  output  stetions  where  time  series  predictions  of  sea  levels 

2.  output  display  options  consisting  of  sea  level  and  current 
plots. 


and  currents  are  stored, 
spatial  field  and  time  series 


The  interfacing  and  operational  environment  for  PCTides, 
between  the  components  of  the  system  and  associated  files,  is 


which  demonstrates  the  relationship 
illustrated  in  Figure  4.3-1 . 


PSI  Technical  Report  SSC-005-00 


PCTides  SDD 


GCOM  running  environment 


Figure  4.3-1 :  Flow  diagram  for  the  GCOM  system. 
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5.0  CSCI  DETAILED  DESIGN 

All  rautines  are  written  in  FORTRAN  90.  The  CSCI  is  platfonn  independent  and  maybe  run  in  the 

^  graphieal  user  interface  (GUI)  interactive^menu  It  can 
also  be  operated  through  use  of  a  command  prompt  in  either  UNIX  or  DOS  mode. 

The  folIo%^mg  paragraphs  give  a  detailed  description  of  the  purpose,  variables,  logic,  and  constraints 
or  tl^  software  el^ente  m  the  CSCI.  Description  of  the  components  of  the  program  will  be 
divided  mto  design  descnptions  of  the  MAPS,  GCOM2D,  and  GCOM3D  models.  ^ 


5.1 


CSC  Maps 


MAPS  IS  a  ^d-pomt,  hydrostalic  primitive  equation  model  capable  of  being  nm  on  varying 

s*"’’®-  Th'  equations  of  motion  are  coded  in  advective  fom 
rather  than  flux  fotm  and  solved  using  a  semi-implicit  time  differencing  scheme  ensuring  that  the 
model  ,s  both  robust  and  eeonomteal  to  run.  even  in  regions  of  steep  terrain.  An  topoZt 
func  ton  of  this  model  ,s  to  provide  high-resolution  representation  of  anemometer  level  S 
and  surface  pressure  as  atmosphenc  boundary  conditions  for  GCOM2D  and  GCOM3D  To  this 
end,  the  turbulence  closure  scheme  has  been  designed  to  allow  the  model  to  be  run  with  its 
ow^t  rnodel  level  at  anemometer  height,  thus  providing  a  direct  simulation  of  the  winds  at  this 
level.  CSC  Maps  consists  of  23  subroutines  and  6  common  blocks,  which  are  interconnected  to 
m^eir  ^  t^lotion  spatial  representation  of  the  wind  and  pressure  field  for  the  GCOM  ocean 


5.7.7 


Constraints  and  Limitations 


5.1.2 


Logic  and  Basic  Equations 


hydrostatic  primitive  equation  model  of  Leslie  et  al  (1985)  It 

thfSix  H  -  solved’rather  than 

the  flux  form.  This  is  done  to  ensure  greater  computational  stability  of  the  model  over  regions  of 

nZ  7""  r  compromising  model  efficiency.  The  governing  equations  for  moiLtum 

Tnica!  ar™  moisture,  utilizing  a  nonnalized  pressure  or  sigma  coordinate  in  the 
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du 

dt 


■  =  -m 


dx  dy 


d(J 


d(p 

ox  dx 


+  (5.1.1) 


dt 


=  -m 


..dv  .dv 

U  — — \-V  -—  -cr  - fU  -  m 

dx  I 


do 


+  RT 


dR 


dtp 

dy  dy 


dtp  ^  RgT 
do 


O 


dW 

do 

dt 


dR 

dt 


=  m 


du  ^ 
ax  a^" 


=  —m 


dx  dy 


- +  —  W -o—W^  =-m  \  U~  +  V- 


O  do  j  do 


+  ant-T\ 


dx  dy  J 


+  Fy+D^,  (5.1.2) 

(5.1.3) 

(5.1.4) 

(5.1.5) 


-m 


3x  dy 


(5.1.6) 


where  U  and  V  are  the  horizontal  wind  components  in  the  direction  of  the  map  co-ordinates  and 

are  defined  to  be  U  =  —  and  V  =  —  where  m  is  the  map  factor,  P,  =  In  p,  where  p.  is  the 

mm 

surface  pressure  and  IV  =  o  +  oW^  where  o  =  — ,  ft  is  the  vertical  velocity  in  o  -coordinates 

P* 

dl\ 

and  .  The  tcnns  tp,  T  and  R  are  the  geopotential  height,  temperature  and  mixing  ratio 

respectively;  and  Fn  are  the  turbulent  exchange  terms  for  momentum,  heat  and 

moisture,  f  is  the  Coriolis  parameter,  ^  is  a  moisture  source  term  and  FI  represents  adiabatic 
heating  through  radiation  and  latent  heat  release;  and  Dr  are  the  horizontal  diffusion 

tenns  for  momentum,  heat  and  moisture  and  Rq  and  Cp  are  the  gas  constant  and  heat  capacity  at 
constant  pressure  for  dry  air. 

5.1.2. 1  Solution  Procedure 


The  equations  shown  above  arc  solved  using  a  semi-implicit  time-diffcrencing  scheme  on  a 
.staggered  .-\rakawa  C-grid  (Messinger  and  Arakawa,  1976)  as  shown  in  Figure  5.1-1.  In  the 
vertical,  the  terrain  following  a -coordinates  system  is  used  as  illustrated  in  Figure  5.1-2,  with  a 
top-down  indexing  system.  The  primary  prognostic  variables,  U.  F.  T  and  R  arc  defined  on  the 
Muir  model  levels  while  a  is  defined  at  the  ‘half  level. 
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Figure  5.1-1:  The  horizontal  grid  structure  used  in  MAPS. 


I’igurc  5.1-2: 1  he  vertical  grid  structure  used  in  MAPS 
full  model  level  except  for  &  . 


■  All  model  \  anables  are  derined  on  the 
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By  introducing  a  reference  vertical  temperature  profile  To(<7)  and  defining  averaging  and 
differencing  operators; 


A'  =\[A{t  +  \M)  +  A{t-\ht)\, 
A’  ^[A{t  +  \M)-A{t-\M)\l 


(5.1.7) 

(5.1.8) 


equations  (5 . 1 . 1 )  -  (5. 1 .6)  can  be  rewritten  using  semi-implicit  differencing  as 

U^'  +  A/|-((r'  -(p,)  +  RcT„  [P^'  A0])= 


where 


=Uit-At)-m^At\  U  —  +V 


dU  . dU 


dx  dy  j  da 


()P  F  D 
ox  m  m 


(5.1.9) 


(5.1.10) 


where 
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R^TAa, 


Bx  B^ 


+  m  - + -  =-mAU — -+V — ^ 


dx  By 


T  W  -  Am  , 

cr  da  \  da 


(5.1.11) 


(5.1.12) 


(5.1.13) 


(5.1.14) 


(5.1.15) 


where 


1.^ 
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'^RHS  ~  A/)  —  m  ^  A/ 
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+ 

+  A/, 


oT  dT 
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f  ccT„  dr 


and 


y  O’  do  J 

dT  dT 

(^~W. -a-^fV. 
do  do 


do 


+  At(H  +  Fj.+Dj.), 


R"  =R(t-At)-m^At\ 


r 


jr'^R  .  dR  ^ 


(5.1.17) 


y 


From  these  equations,  the  hydrostatic  equation  is  solved  to  provide  geopotential  heights  The 
thermodynamic  equation  is  then  transformed  into  a  Helmholtz  equation  using  the  contln^itv  and 

s™™  “"'“‘y  “0  equates”" 


5.1.2.2 


Initialization 


The  model  is  initialized  by  nesting  within  interpolated  analysis  and  prognosis  fields  from  a 

coarser  gnd  model  covenng  a  large  region  (usually  the  globe).  In  this  case  the  model  is  nested 

Atmospheric  Prediction  System  (NOGAPS)  (Hogan  et 
al.,  1991)  available  at  0000,  0600,  1200  and  1800  UTC  each  day.  ^  ^  ^ 


5. 1.2.3 


Boundary  Conditions 


The  bound^y  conditions  are  applied  using  a  nesting  procedure  that  blends  the  extemallv 
prescnbed  fields  with  the  model-simulated  values  over  a  number  of  grid-points  the 
computational  domain.  Absolute  boundary  conditions  are  applied  to  each  model  variable  7 
with  decreasing  intensity  from  the  boundary  to  some  specified  number  of  grid-points  (typically 
less  than  or  equal  to  10)  within  the  domain,  according  to  the  following  equate 


d  -  (I  -OC)<pp  3-Ot(P^ 


(5.1.18) 


where  ,s  the  model  predicted  x  alue,  is  the  externally  applied  %  alue,  provided  by  analyses 

in  this  case,  and  a  varies  from  1  at  the  boundary  to  zero,  n„ 
to  a  cosine  function  ’  " 


points  into  the  domain  according 


a  =  0.5(cos  ;r ( 1  -  /;  / =  1  „ 


(5.1.19) 


The  penultimate  grid-points  are  then  smoothed  to  dampen  out  graMt>'  wa%  es  that  will  otherwise 
accumulate  at  the  boundaries  and  lead  to  numerical  instability. 
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The  model  is  capable  of  being  nested  within  itself  to  provide  higher  resolution  simulations. 
Under  these  circumstances,  the  boundary  conditions  are  applied  the  same  way  as  described 
above  where  the  are  provided  by  a  previously  performed  lower  resolution  run. 

5. 1.2.4  Temporal  Filtering 

A  temporal  filter  is  applied  to  all  prognostic  variables  to  reduce  the  amount  of  energy  in  waves 
with  high  frequencies.  All  variables  are  filtered  according  to 


2  ^  (5.1.20) 

(Asselin,  1972)  where  ^  is  the  filtered  variable.  The  coefficient  i/has  a  value  of  0.8. 

5. 1. 2. 5  Parameterization  of  Sub-grid  Scale  Processes 

Terms  go\  eming  the  effects  of  various  processes  such  as  turbulent  vertical  diffusion,  radiation 
and  precipitation  are  included  in  the  prognostic  equations  for  the  horizontal  wind  components, 
temperature  and  moisture.  Some  of  these  processes  can  be  calculated  explicitly  while  others 
must  be  parameterized.  In  the  model,  the  effects  of  these  various  processes  are  treated  as  a 
separate  adjustment  to  the  adiabatic  semi-implicit  solution  of  the  equations  of  motion. 

5. 1.2.6  Turbulence  Closure  Scheme 

This  section  describes  the  representation  of  the  vertical  transports  via  subgrid-scale  eddies.  These 
turbulent  flux  tenns  arise  as  dissipative  terms  in  the  equations  of  motion  through  Reynolds 
averaging  considerations.  In  the  V  and  F-momentum  equations,  the  significant  contribution 

comes  from  the  vertical  components  of  the  stress  tensor,  -  puw  and  -  pvw,  respectively,  while 

in  the  thermodynamic  and  moisture  equations  the  analogous  terms  -  pOw  and  -  pWv  represent 
the  vertical  transfer  of  heat  and  moisture,  respectively.  These  flux  terms  give  rise  to  a  closure 
problem  that  is  generally  overcome  by  parameterizing  in  terms  of  mean  flow  variables. 

The  turbulence  closure  scheme  available  in  MAPS  is  a  simplified  second  order  closure  scheme 
that  is  described  by  Mclnnes  and  Hess  (1992).  By  definition,  second  order  closure  implies  that 
the  model  carries  prognostic  equations  for  seeond  order  moments  and  closure  of  the  system  of 
equations  is  achieved  by  parameterizing  or  simplifying  third  order  moment  terms.  However, 
Mellor  and  Yamada  (1974)  describe  a  hierarchy  of  second  order  closure  models  of  varying  levels 
of  complexity  ranging  from  the  most  complex  “level  4”  model  to  the  least  complex  “level  2” 
model.  In  simplest  tenns,  the  higher  the  level,  the  greater  the  number  of  prognostic  equations  for 
the  second  order  moments  that  arc  retained  in  the  scheme. 

In  MAPS,  the  turbulence  closure  scheme  is  the  so-called  “level  2'A"  model  described  in  Galperin 
et  al.,  (1988),  which  is  similar  in  complexity  to  the  “level  I'A"  model  of  Mellor  and  Yamada 
(1982).  In  this  model,  a  prognostic  equation  is  retained  for  the  turbulent  kinetic  energy,  q,  only, 
while  all  other  turbulent  flux  equations  arc  subjected  to  various  simplifications  and 
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approximations  to  reduce  them  to  a  set  of  diagnostic  equations 
and  mean  flow  variables. 


that  relate  other  turbulent  fluxes 
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(5.1.21) 

(5.1.22) 

(5.1.23) 

(5.1.24) 

(5.1.25) 

(5.1.26) 

(5.1.27) 

(5.1.28) 

(5.1.29) 

(5.1.30) 


The  equations  anse  from  the  Reynolds  Stress  equations  after  scaling  assumptions  have  been 
made  and  par^unetenzat.ons  introduced  for  many  of  the  tenns  such  as  the  triple  Lrrelation  terms 


( A, ,  A ,  5, ,  ,  C, )  =  (.92,. 74, 1 6.6,1 0. 1,.08). 


(5.1.31) 


fhc  turbulent  length  scale  /  which  forms  the  basis  of  the  second  order 
detenmned  using  an  algebraic  expression  suggested  by  Blaekadar  (1962) 


closure  assumptions  is 
abo\'e  the  surface  layer; 


/  = 


\ 

y 


(5.1.32) 
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where  k  is  von  Karman’s  constant  (/:=0.4)  and  z  is  the  height.  Mellor  and  Yamada  (1974)  specify 
X  according  to; 


A  = 


(5.1.33) 


jqdz 

0 

where  a  is  an  empirical  constant. 

These  equations  are  used  to  solve  for  the  vertical  eddy  flux  terms  -  and  which  in  turn 
allow  eddy  coefficients  to  be  specified  using; 


—  UW  —  vw 


U- J 

[dzj 

^  ~0w 
i  9z 


(5.1.34) 


(5.1.35) 


These  eddy  coefficients  are  applicable  at  the  center  of  each  model  layer.  The  consequent 
subgrid-scale  source  terms  in  the  momentum  and  thermodynamic  equations  in  the  model  are  then 
simply  obtained  by  the  vertical  divergence  of  these  fluxes  in  the  \  ertical  diffusion  equations 
which  for  the  u-component  of  velocity  takes  the  form; 

dr  dz  9z  (5.1.36) 

The  second  order  closure  scheme  is  applied  at  all  levels  above  the  surface  layer.  The  eddy 
coefficients  that  arise  from  the  algebraic  manipulation  of  equations  (5.1 .21-5.1.30)  are; 


IqAi 


1-3C,  -^-3A^G„ 


B 


(52-3T,) 


1- 


6A, 


B 


lJ 


-3C^{6A^  +B2) 


(l-9A,A,Gf,  )il-3A,G„i6A,+B,)) 


B 


■  M 


U. 


(5.1.37) 


(5.1.38) 


1  -  3A2G n  (6/l|  +  Bj ) 

,  „  (I  q  330 

where  G,  =  -  —  —  — 

■■  U/  ^  dz  ' 


1  he  prognostic  equation  for  the  turbulent  kinetic  energy  is; 
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~du  —dv  g  - 

-uw— - vw—-¥  —  w6^-~, 

dz  0,  B,/ 


(5.1.39) 


where  the  constants  and  B3  are  0.2  and  0.3  respectively,  -  sj  ^ 

V  P  ^ 

gradient  heat  flux,  and  is  the  convective  velocity  defined  as  follows- 


(-h\^  . 

n;  =  J  “*  —  I  ,  in  convective  conditions; 
otherwise. 


'eg  -  ^  “t; - r  is  the  counter 

\pCpW.h\ 


(5.1.40) 


ertical  diffusion  equation  is  solved  for  the  horizontal  wind  components,  temperature  and 
mixing  ratio  using  an  implicit  numencal  finite  difference  scheme  which  is  given  below  for  U- 

L"'*' -U''  ..  rt/r,‘ +1/-' ^ 


=  K  ^^k-\ 

""  2Az 


(5.1.41) 


5.1.2.7 


Monin-Obukhov  Surface  Layer 


'T  l'  “.f"®  similarity  theory  (Monin  and 

Obukhov  1954  .  This  layer,  which  is  typically  about  10  m  deep,  is  assumed  to  lie  between  the 

.h: 


dlV  u,  fz) 

3--  fe  \l} 

fz^ 

A’pW.fe  "  L  \' 


dR  E  { z 
—  =  — q)  _ 

3-  pa.kz  "\l 


(5.1.42) 


(5.1.43) 


(5.1.44) 


where  ff 


-)  2  — 

(f/‘  +  F  )" ,  Cp  is  the  specific  heat  of  dry  air. 


the  air  density  and  k  is  the  Monin-Obukhov  length  defined  as- 


is  the  roughness  length,  p  is 


(5.1.45) 


where  0,  is  the  virtual  potential  temperature,  //  is  the  sensible  heat  flux.  The  eddy  viscosity 
cocfilcients  arc  defined  accordine  to; 
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K„  =M./0„  -  . 


(5.1.46) 

(5.1.47) 


According  to  Dyer  (1974),  the  following  stability  functions  are  found  to  provide  the 


accurate  flux -gradient  description;  in  the  unstable  case  [  -^  ]<  0 

UJ  ’ 

=fl-16f-ll  \ 


z  \\  ^ 
■16  -  , 

U 


(5.1.48) 

(5.1.49) 


and  for  the  stable  case  —  >  0 

L 


=  1+5 


(5.1.50) 


In  neutral  conditions,  -  ‘I’ at  -  1  giving  rise  to  the  well  known  logarithmic  wind  profile  ir 

equation  (^.  1 .42).  Equations  (5. 1 .48)  and  (5. 1 .49)  can  be  integrated  analytically  to  give; 


i  -  =V,  =—  In 

K  max  ^  k  max  ,  i 


j  ^max  ,  VT/  ^0 

^\~irv^^AT 


(5.1.51) 


©in.x-0..=  In 

pCpU.k  I 


where 


~  ^  _  VP  I  ^kuMX  ,  wy  f 

~  — 

-T  L 


21n  kin  -2tan-'cl>;;,  ^  <  0; 


'i^tr(0=  0, 

-5^, 


^>0. 


(5.1.52) 


(5.1.53) 


(  1  -f  (J)  '  A 
21n  "  . 

2 

,  C<0; 

(5.1.54) 

0, 

r  = 

V 

O 
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The  roughn^ess  length  zo,  over  the  land  is  assumed  (following  eommon  practice)  to  be  a  physical 

foimulLiL  of 

Chamock  (1955)  where  the  roughness  length  is  a  function  of  the  wind  speed- 
IrJ 

(5.1.55) 


Zq  =  .032^ 


t^r/=  It temperature  equation,  zr,  is  related  to  the  roughness  length  according 
Surface  Temperature 


to  Zj-—  .1  Zq. 

5. 1.2.8 


Surface  ternperature  over  land  masses  is  calculated  using  a  3-layer  slab  model  based  on  the 

Medium  Range  Weather  Forecasting  (ECMWF)  model, 
(ECM\W,  1984).  The  model  is  illustrated  in  Fig.  5.1-3  where  Ts  is  the  surface  soil  temperature 
assumed  to  v^y  on  a  diurnal  time  scale.  To  is  the  deep  soil  layer  temperature  assumed  to  vary  on 
a  time  scale  of  several  days.  Ta  is  assumed  to  be  sufficiently  deep  that  the  temperature  varies  on 

climatological  time  scales  and  as  such  can  be  held  constant.  The  depths  of  the  layers  are  defined 
in  cm  as  follows; 


(D,,£)2,Z)3)  =  (7.2,43.2,43.2). 


(5.1.56) 


Di 


Rs  Rl  h  le 

1 1 1 1 

{ 


D2< 


D3 


I 


-Ta- 


igure  5.1-3:  Schematic  diagram  showing  the  3-laycr  surface  temperarure  prediction  scheme 
1  he  quantities  Rs.  R,_.  H  and  LE  represent  the  incoming  shortwave  and  longwave  radiation 
sen.sible  and  l^atcnt  heat  lluxes,  respectively.  /Ivand  To  are  predicted  mean'^temperatures  oCer 
depths  D,  and  D.  respectively  and  Ta  is  held  constant.  The  thickness  of  the  depths  D,.  D,  and 
D?  arc  7.2,  43.2,  and  43,2  cm,  respectively. 
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The  equations  for  the  surface  and  deep  soil  layers  are; 
dT,  _  IF  ^ 

dt  \D,{D,+D^y  (5.1.57) 

dT,  _  {T,-T,)k  ^(Ta-T,)K 

dt  iD,(D, +£>2)  ’  (5.1.58) 


where  is  the  soil  density,  Q  is  the  specific  heat  capacity  of  the  soil  and  k  is  the  thermal 

diffusiviU'.  The  transfer  of  heat  through  the  soil  is  assumed  to  be  through  diffusive  processes 
except  in  the  surface  layer  where  the  first  term  on  the  RHS  of  equation  (5.1.57)  represents 
radiative  and  convective  transfer  of  temperature  with  the  atmosphere: 

1F  =  R,+R,+Es+H, 

(5.1.59) 

and  Rs  is  the  short  wave  solar  radiation,  Rl  is  the  long  wave  radiation  and  £5  and  H  are  the  latent 
and  sensible  heat  respectively  defined  according  to; 


//  =  pc^c,,|u*_|(r,_-7;), 
and 

^ s  ~  P^E  ^ n  ^  |(?jo,  ~  Ik  niax  )> 


(5.1.60) 

(5.1.61) 


where  Ch  is  the  heat  transfer  coefficient,  is  the  latent  heat  and  a  is  an  efficiency  coefficient 
defined  as  ,  where  W,  is  the  soil  moisture  and  is  the  value  of  when 

saturation  occurs.  A  simplified  surface  radiation  scheme  is  used  to  incorporate  the  effects  of  the 
radiative  flux  on  surface  temperature  variation.  This  scheme  provides  the  shortwave  radiation  in 
cloudy  and  clear  sky  conditions  and  a  long  wave  radiation  budget  for  frequencies  between  800 
and  1200  s’  .  Complete  cooling  to  space  is  assumed  for  all  frequencies  outside  this  range. 

Equation  5. 1 .57  is  solved  using  the  following  finite  difference  procedure; 


-J-  n  +  \  _  /J 


At 


(5.1.62) 


where  0(7^)  repiesenls  the  right  hand  side  of  equation  5.1.57.  B\’  expanding  this  term  in  a 
Taylor  senes  about  <p{TX  this  equation  becomes 
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which  simplifies  to  the  following  explicit  relationship  for  ; 


+  I  _  rnri  _ 

^  ~T 


1-A^ 


V 


dT 


n  \ 


(5.1.64) 


This  ‘Euler  Backward’  procedure  is  required  only  for  since  once  this  is  known,  T, 
solved  explicitly. 


d  can  be 


5.L2.9 


The  Hydrological  Cycle 


A  3-layer  sub-surface  model  (Figure  5.1-4)  analogous  to  the  scheme  used  for  surface  temperature 
prediction  is  used  for  the  calculation  of  surface  wetness.  The  governing  equations  are; 


Ph,o  "  \D,{D,+D^y  (5.1.65) 

dW,  _  (W^-W,)X  ^(W^,-W,)X 

dt  jD^iD^+D^)  ’  (5.1.66) 


where  E  and  Pr  represent  the  moisture  transfer  to  the  ground  from  water  vapor  and  precipitation 
respectively  and  and  A  are  the  density  and  diffusivity  of  water.  The  values  of  Z),,  Dj  and 

Di  are  the  same  as  those  used  for  the  surface  temperature  prediction  scheme  and  the  moisture  in 
the  lowest  sub-surface  layer  fVa  is  specified  using  climatological  values. 

Pr  E 


Di  ^  - - - Runoff 


I 


Figure  5.1-4:  Schematic  showing  the  3-laycr  moisture  scheme  used  by  the  model.  The 
quantities  E  and  PR  represent  the  evaporation  and  precipitation  respecti\  ely,  while,  W,  and  Wj 
are  the  predicted  moisture  content  over  the  depths  D|  and  Do.  * 
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5. 1.2.10  Large  Scale  Precipitation 

Large-scale  condensation  occurs  in  the  model  when  the  mixing  ratio  R  exceeds  a  specified 
fraction  of  the  saturated  mixing  ratio  Rsat-  The  condensation  rate  is  gi\'en  by 


2M{\^-L\^Jc^RJ*^y  (5.1.67) 

Where  T*  and  i?*  are  the  predicted  values  of  temperature  and  mixing  ratio  before  condensation. 
Due  to  the  fact  that  latent  heat  release  during  condensation  increases  both  temperature  and  Rsat, 
not  all  of  the  excess  water  condenses.  The  final  values  of  R  and  T  are  given  by; 

R  =  R*-2StC  (5.1.68) 

T  =  T*+~2MC 

c,  (5.1.69) 


Condensed  water  is  assumed  to  fall  out  immediately  with  no  subsequent  evaporation  occurring. 
The  precipitation  rate  due  to  large-scale  rain  is  then  given  by; 


P<  =  — 


1 


k  ntax 


Y,CiP*Aaj, 


(5.1.70) 


where  C/,  is  the  condensation  rate  at  level  k. 


5.1.2.11  Cumulus  Convection 

Subgrid-scalc  cumulus  convection  is  parameterized  using  a  scheme  based  on  Kuo  (1965)  with 
modifications  suggested  by  Hammarstrand  (1977).  Cumulus  convection  is  assumed  to  occur 
when  low  level  moisture  convergence  occurs  in  conjunction  with  an  upper  layer  of  conditionally 
unstable  air.  The  extent  of  the  cloud  is  taken  to  be  the  region  bounded  by  the  lifting 
condensation  level  of  air  near  the  surface  and  the  intersection  of  the  moist  adiabat  with  the 
environmental  sounding  at  upper  levels.  Cumulus  clouds  formed  in  this  wa\'  arc  assumed  to  mix 
instantaneously  and  completely  with  the  environment.  The  moistening  and  heating  rates  arc 
proportional  to  the  differences  between  the  cloud  temperature,  Tc  and  mixing  ratio,  Rc  and  the 
environmental  values  7’ and  R. 


The  cumulus  convection  scheme  is  divided  into  .several  steps; 

(i)  Check  the  sign  of  the  large-scale  moisture  convergence,  where 
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M 


CO/IV 


^UR  dVR  dcR 

-^  +  -z—  + - 

ox  dy  dz 


\ 

da 


If  Afconv>0  then  proceed  to  next  step. 

(ii)  Check  for  conditional  instability  in  the  model  vertical  profile.  That  is,  if 


(5.1.71) 


a  , 

-—{(p  +  CpT  +  LR)<0 


(5.1.72) 

(5.1.73) 


(iii)  Calculate  cloud  base  by  carrying  out  a  dry  adiabatic 
model  level. 

(iv)  Calculate  the  cloud  temperature  profile  Tc(a)  from  the 

ac 


da 


ascent  of  air  starting  from  the  lowest 
equation  for  the  moist  adiabat; 

(5.1.74) 


fs  constant  for  water  vapor.  This  equation  is  integrated  using  a  Runge-Kutta 
method  and  is  halted  when  cloud  top  is  reached.  ^  ^ 

(v)  The  heating  and  moistening  at  level  k  are  assumed  to  be  proportional  to  the  differences 
between  the  model  cloud  and  the  environment, 

^T,=^riTc-T) 


=^ARc-R) 


(5.1.75) 


It  is  assumed  that  ^^.-^and  where/? is  a  positive  constant  with 
g  is  calculated  from 


(5.1.76) 

value  of  0. 1 .  The  value  of 


where 

=~  ]cAT,-T)da 

^  CT/ 

p  *  7 

Q:  = -  \L(R  R),JC7 

Finally,  the  convective  precipitation  is  computed  from 
p  *  7  c 

Pc  =pg —  \  —(l\.-T)da 
.‘t  i  i- 


(5.1.77) 

(5.1.78) 

(5.1.79) 


(5.1.80) 
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5. 7.2, 12  Shallow  Convection 

A  simple  shallow  convection  scheme  is  employed  up  to  level  o=0.1.  This  scheme  was  added  to 
counter  the  tendency  of  the  deep  convection  scheme  to  cause  too  much  subsidence  and  moisture 
accumulation  in  the  lowest  model  layers,  particularly  in  the  tropics.  Shallow  convection  is 
assumed  to  occur  if  instability  is  observed  in  the  levels  above  or  below  and  the  relative  humidity 
of  the  current  level  is  greater  than  75%.  In  the  event  of  shallow  con\'ection  occurring,  the  eddy 
viscosity  coefficients  of  momentum  are  increased  to  reflect  the  greater  turbulence  occurring 
within  the  level. 

5.1.3  Data  Elements 

5.1. 3.1  Commonly  Used  Parameters 


File;  amparam.for 


Variable 

Type 

Description 

nk 

Integer 

Maximum  vertical  dimension 

ni 

Integer 

Maximum  horizontal  dimension  (x-direction) 

jii _ 

'  Integer 

Maximum  horizontal  dimension  (v-direction) 

jiii _ 

^  Integer 

NI*NJ 

These  parameters  are  used  in  subroutines  AMOUT,  ASSLIN,  BNDFIX2  CONST  CYCLE 
FILTER,  FILTER,  FIZZIX,  GETDT,  PUTOUT,  INNER,  INNER2,  NESTIN  PARAMs’ 
SEMIMP. 


File:  amcommon.for 


Variable  i  Type 

Description 

nkp  1  Integer 

NK+l  '  ~ 

nout  integer 

Number  of  output  stations 
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5.1.3.2  Internally  Defined  Parameters 


The  following  parameters  are  used  in  subroutine  FIZZEX. 


Variable 

Type 

Description 

yal,ya2,ybl,  yb2,vb3.vcl 

Real 

Constants  used  in  turbulence  closure  scheme 

zal,za2,  za3,  za4 

Real 

Constants  used  m  solution  of  surface  temperature  and' 
moisture  equation. 

yml 

Real 

Empirical  constant,  cl,  used  in  calculation  of  mixing  length 
in  turbulent  diffusion  scheme. 

el,  e2,  e3 

Real 

Constants  used  m  solution  of  vertical  diffusion  equation 

cq 

Real 

LConstant  used  m  the  calculation  of  (9  .„fi-om  6  and/? 

emissg,  emissc,  stefbo 

Real 

Ground  emissmty,  cloud  emissmty,  and  Stephan-Boltzman 
constant. 

zrl,  ztl 

Real 

Constants  used  in  calculation  of  cloud  transmissiviry 

blimit,  rhkuo 

Real 

Highest  allowed  cumulus  base,  critical  relative  humidity 
(Kuo  scheme). 

rv,  hi 

Real 

Constants  used  in  large-scale  rain  scheme 

The  following  parameters  are  used  in  subroutine  NESTIN. 


Variable 

Type  1 

- -  V  A  X  Xi  >(  . 

Description 

nbr,  nbrp 

Integer  | 

Dimension  for  nestmg  fields,  NBRP=NBR  +  1 

Subroutines 


5.1. 3.3  Subroutine  AMOUT 


The  purpose  of  this  subroutine  is  to  output  model  fields  in  direct  access  fonnat 
MAIN  program. 


It  is  called  by  the 


Calling  sequence:  subroutine  amout(lu2,lu5,dt) 

Data  declarations:  real  dt 

logical  Iu2,lu5 

Arguments:  lu2  Logical  unit  number  of  the  model  run  log  file. 

Logical  unit  number  of  the  direct  access  file. 
Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 
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5.L3.4  Subroutine  ASSLIN 

ASSLIN  performs  temporal  filtering  of  U,  V,  T  and  RM  arrays  and  is  called  by  the  main  program. 
Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

Algorithm; 

Refer  to  Temporal  Filtering,  section  5. 1.2.4. 

5.L3.5  Subroutine  BNDFIX2 

BNDFIX2  smoothes  penultimate  grid  values  of  arrays  PSP,  UP,  VP,  RMP,  TP,  OMEGA,  TSP, 
WSP,  TDP  and  WDP  after  lateral  boundary  conditions  are  applied  to  prevent  accumulation  of 
gravity  wa\'es.  Called  by  the  MAIN  program. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

Algorithm: 

Refer  to  5. 1 .2.3,  Boundary  Conditions. 

5. 1.3.6  Subroutine  CONST 

CONST  sets  up  constants  for  horizontal  diftusion,  and  eigenvectors  and  matrices  for  semi-implicit 
solution  procedure.  It  is  called  by  the  MAIN  program. 

Calling  sequences;  subroutine  const(dt) 

Data  Declarations;  real  dt 

Arguments:  dt  Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5. 1. 3. 7  Subroutine  CYCLE 

CYCLE  rotates  time-stepping  arrays  at  the  end  of  each  time  step.  It  is  called  by  the  MAIN  program. 
Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5. 1.3.8  Subroutine  DADADJ 

Subroutine  DADADJ  perfonns  dry  adiabatic  adjustment.  It  is  called  by  SEMIMP. 
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*. 


Calling  sequences:  subroutine  dadadj(dq,  dqmqa5,tx,  nz,  ipc,  jpc) 

Data  Declarations:  real  dq,  dqmqaS,  tx 
integer  nz,  ipc,  jpc 

Arguments:  dq  array  of  A  a. 

dqmqaS  Convective  adjustment  array, 

be  Temperature  calculated  by  DAD ADJ. 

Vertical  dimension. 

jP*'  I-grid  point  (for  diagnostic  print). 

jP^  j-grid  point  (for  diagnostic  print). 

5J.3.9  Subroutine  FILTER 


performs  spatial  filtenng  of  multi-level  variables  or  their  tendencies.  Single  level  fields  are 
filtered  after  eveiy  physics  time  step  NPHYS.  It  is  called  by  the  MAIb,^  program. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 
5.1.3.10  Subroutine  FIZZIX 

FIZZIX  calculates  adjustments  to  U,  V,  T  and  RM  due  to  vertical  diffusion,  cumulus  convection 

“calwt; 


Calling  sequences:  subroutine  fizzix(lu2,  luO,  dt) 

Data  Declarations:  real  dt 

logical  lu2, luO 

Arguments:  lu2  Logical  unit  number  of  the  model  run  log  file. 

foO  Logical  unit  number  of  the  direct  access  file. 

Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 
5.1.3.  / /  Subroutine  GETDT 


GETDT  e\  aluatcs  the  time-step  based  on  tlic  CFL  condition  and  is  called  by  the  MAIN  program. 
Calling  sequences:  subroutine  getdt(lu2,  lu3,  ntime,  spmax,  dt) 
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Data  Declarations:  real  dt,  spmax 
integer  ntime 
logical  lu2,  lu3 


Arguments:  lu2 
lu3 
ntime 

spmax 

dt 


Logical  unit  number  of  the  model  run  log  file. 
Logical  unit  number  of  the  direct  access  file. 
Counter  incremented  by  1  each  time  a  nesting  file 
exists. 

Maximum  wind  speed  found  in  the  wind  field  array. 
Model  time  step. 


Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 


5.13.12  Subroutine  PUTOUT 


PUTOUT  outputs  mn-time  information  to  the  screen  either  in  text  fomiat  or  graphically.  It  also 

wntes  a  run-time  log  file  (MAPS.LOG)  and  the  stations  time  series  file.  It  is  called  by  the  MAIN 
program. 

Calling  sequences:  subroutine  putout(llog,  lutsd,  lutsp,  lu2,dt) 

Data  Declarations:  real  dt 

logical  Hog,  lutsd,  lutsp,  lu2 


Arguments:  Hog 

Logical  unit  ofMAPS.LOG. 

lutsd 

Logical  unit  of  the  TSD  files. 

lutsp 

Inacti\'e  parameter. 

lu2 

Inactive  parameter. 

dt 

Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 


5. 1.3.13  Subroutine  INNER 


INNER  calculates  the  gcopotential,  calls  INNER2,  and  calculates  the  forcing  tenn  for  the  solution 
of  the  Helmholtz  equation.  It  is  called  by  the  MAIN  program. 

Calling  sequences:  subroutine  inner(dt) 

Data  Declarations:  real  dt 

Arguments:  dt  Model  time  .step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 
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5.1.3.14  Subroutine  INNER2 

INNER2  calculates  the  right  hand  side  of  the 
called  by  subroutine  INNER. 


semi-implieit  equations  for  U,  V,  T,  and  RM. 


It  is 


Calling  sequences:  subroutine  inner2(dt) 

Data  Declarations:  real  dt 

Arguments:  dt  Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5. 1.3. 15  Subroutine  NESTIN 


NESTIN  on  imtial  call  reads  direct  access  header  infomiation,  sets  sigma  levels  and  TEAR 
zmhahzes  gnd  projechon  ^d  calculates  grid  points  nearest  to  specified  stations  for  time  se^s 
outpu  T^e  weights  are  calculated  for  the  nesting  scheme,  as  well  as  solar  declination  and  the  top 
yer  at  which  vertical  diffiision  is  applied.  Convective  adjustment  variables,  constants  for  the  semi- 
imphcit  scheme  and  map  factors  at  staggered  grid  points  are  defined.  On  subsequent  calls 
bound^  ^ndihons  are  applied  to  model  predicted  variables  via  a  nesting  tech^que  TWs 
subroutine  is  called  by  the  MAIN  program.  ®  ^ 


Calling  sequences:  subroutine  nestin(lstn.  Hog,  lu2,  lu3,  dt,  iend) 

Data  Declarations:  real  dt 

integer  iend 

logical  Istn,  Hog,  lu2,  lu3 


Argument:  Istn 

Hog 
lu2 
lu3 
dt 

iend 


Logical  unit  number  of  the  model  run  log  file. 
Logical  unit  number  of  the  model  run  log  file. 
Logical  unit  number  of  the  model  run  log  file. 
Logical  unit  number  of  the  direct  access  file. 
Model  time  step. 

Nesting  file  flag.  Set  to  one  after  last  nesting  file  is 
read. 


Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5. 1.3.16  Subroutine  P ARAMS 


This  subroutine  reads  the  MAPS.DAT  file  and  sets  various  constants  and 
the  MAIN  program. 


parameters. 


It  is  called  by 


Calling  sequences:  subroutine  const(lpar,  lu2) 
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Data  Declarations:  logical  Ipar,  lu2 

Arguments:  Ipar  Logical  unit  number  of  the  model  run  log  file. 

Iu2  Logical  unit  number  of  the  direct  access  file. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5. 1.3.17  Subroutin  e  SEMIMP 

SEMDVTP  sets  up  the  Helmholtz  equation  to  solve  for  TP  and  also  solves  for  UP,  VP,  and  PSP.  It  is 
called  by  the  MAIN  program. 

Calling  sequences:  subroutine  semimp(dt) 

Data  Declarations:  real  dt 

Arguments:  dt  Model  time  step. 

Common  Blocks:  DOUBLE,  THREED,  TWOD,  ONEDR,  ONEDC 

5.1.3.18  Subroutine  TRIDAG 

TRIDAG  is  a  tridiagonal  matrix  solver.  It  is  called  by  FIZZIX  to  solve  the  vertical  diffusion 
equation. 

Calling  sequences:  subroutine  tridag(a,  n,  ndim) 

Data  Declarations:  real  a 

integer  n,  ndim 

Arguments:  a  a(ndim,  4)  is  the  array  to  be  inverted. 

n  The  number  of  elements  stored  in  first  dimension  of 

a. 

ndims  The  dimension  of  the  first  array  element  of  a. 

5.1.3.19  Function  FF 

FF  is  used  in  the  calculation  of  tlie  Monin-Obukhov  length  L.  It  is  called  by  FIZZIX. 

Calling  sequences:  function  ff(zl,  zOI) 

Data  Declarations:  real  zl,  zOl 
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Arguments:  zl 
zOI 

5.1.3,20  Function  GG 


Lowest  layer  depth  divided  by  Monin-Obukhov 
length. 

Roughness  length  (zO)  divided  by  Monin-Obukhov 
length. 


GG  IS  used  m  the  calculation  of  the  Monin-Obukhov  length  L.  It  is  called  by  FIZZIX. 
Calling  sequences:  function  gg  (zl,  zhl) 

Data  Declarations;  real  zl,  zhl 


Arguments:  zl 
zhl 

5.1.3.21  Function  PSIM 


Lowest  layer  depth  divided  by  Monin-Obukhov 
length. 

zh  divided  by  Monin-Obukhov  length. 


s^celayr^culata'^  f"  Monin-Obakhov 

Calling  sequences:  function  psim  (zl) 

Data  Declarations:  real  zl 


Arguments:  zl 

5.1.3.22  Function  PSIH 


Lowest  layer  depth  divided  by  Monin-Obukhov 
length. 


This  is  a  profile  stability  function  for  heat, 
calculation. 


It  is  called  by  FIZZIX  for  Monin-Obukhov  surface  layer 


Calling  sequences;  function  psih  (zl) 
Data  Declarations;  real  zl 


Arguments:  zl 

5.1.3.23  Function  PlllMI 


Lowest  layer  depth  divided  by  Monin-Obukho\' 
length. 


PHIMI  calculates  4,,,  for  Monin-Obukhov  calculation.  It  is  called  by  FIZZIX. 
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Calling  sequences:  function  phimi  (zl) 

Data  Declarations:  real  zl 

Arguments :  zl  Lowest  layer  depth  divided  by  Monin-Obukho v 

length. 

Algorithm: 

Refer  to  Monin-Obukhov  Surface  Layer,  Section  5. 1.2.7. 

5.13.24  Function  PHIHl 

PHIMI  calculates  for  Monin-Obukhov  calculation.  It  is  called  by  FLZZIX. 

Calling  sequences:  function  phihi  (zl) 

Data  Declarations:  real  zl 

Arguments:  zl  Lowest  layer  depth  divided  by  Monin-Obukhov 

length. 

Algorithm: 

Refer  to  Monin-Obukhov  Surface  Layer,  Section  5. 1.2.7. 

5.1.3.25  Function  ESAT 

ESAT  calculates  the  saturated  vapor  pressure  for  a  given  temperature.  It  is  called  by  FLZZIX. 
Calling  sequences:  function  esat  (t) 

Data  Declarations:  real  t 

Arguments:  t  Temperature  for  which  the  saturated  vapor  pressure  is 

required. 

The  organization  of  subroutines  and  functions  within  the  MAPS  model  is  illustrated  in  Ficure 
5. 1.3-1.  ^ 
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Main  Program 


Subroutines 

.  PARAMS 

Read  MAPS.DAT  and  set  flags 

NESTIN 

:>et  titnestep,  boundarytendencies 
CYCLE 

Assign  values  to  timestep  arrays 


Subroutines/Functions 


GETDT 


►j _ PUTOUT 

►  AMOUT 

WritcQut  MAPS.OOO  file 

p  CONST 

Set  up  constants,  eigenvectors  etc 


EIGENP 


NtATINV 


Time¬ 

stepping 


_ ADATE2 

INNER 

Solves  semi-implicit  equations 
SEMIMP 

Solves  helmholtz  equations 
FIZZIX 

_ Physical  processes 

NESTIN 

Applies  boundary  conditions 
BNDFIX2 

Smooths  after  nesting 

FILTER 
_ Filter  fields 

ASSLIN 


OPTIMUM 


CYCLE 

Rotates  time-stepping  arrays 
PUTOUT 

Outputs  mn-time  info 
AMOUT 

Writes  to  direct  access  file 


TRIDAG 


D  \  routines 


Figure  5.1  3-1:  Flowchart  illustrating  the  mnning  of  the  MAPS  model.  Unshaded  subprograms 
indicate  those  stored  in  separately  compiled  libraries.  Note  also  that  multiple  calls  to 
subprograms  within  subroutines  are  not  indicated  to  simplify  the  flow  diagram. 
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5. L4  The  MAPS  System 


In  addition  to  the  MAPS  model,  several  other  modules  and  programs  are  necessary  to  complete 
the  modeling  system.  These  include  programs  for  setting  up  model  grids  over  specified 
geographic  regions,  extraction  and  conversion  programs  to  pro\'ide  atmospheric  initial  and 
boundary  conditions  in  suitable  format  for  MAPS  and  display  programs  for  viewing  model 
output.  These  additional  programs  are  summarized  in  Figure  5. 1.4-1  while  the  key  additional 
elements  are  described  in  the  following  sections. 


MAPS  running  environment 


User  Specified 

Inputs  - - 

(GRIDGEN.DAT) 


Grid  Generator 
GRIDGEN 


TOPOG.DAT 


User  Specified 

Inputs  - 

(MAPS.DAT) 


User  Specified 
Inputs  _ 


(M  A  PS. DAT) 
(STATIONS.DAT) 


User  Specified 

Inputs  - 

(via  menu) 


Archived  Topographic  Data 


(NESTING) 


NOGAPs.ooo  User  Specified  Inputs 


_ 3 

^ - X - NOG  APS.006  CIC _ fM  APS  D  AT) 

Field  G< 

AMSU 

i 

snerator 

BGRD 

◄ - 

Input  Conversion  Program 

NOGAPS 

AM  SUBG  RD.OOO 
AMSUBGRD.006  etc 


Preprocessor 

AMPREPSL 

AMPREPPL 


Archived  Climatological 


Data  I 


AM  PREP.OOO 
AM  PREP.006  etc 


Mesoscale  Model 
MAPS 


M  APS.OOO 
M  APS.006  etc 


(NESTING) 


STATIONl.TSD 
STAT10N2.TSD  ci 


Display  Module 
AMD  ISP 


Display  Module 
TSPLOT 


Figure  5.1.4-1:  Flowchart  illustrating  the  running  of  the  .MAPS  system. 
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5J.4.1 


The  Grid  Generator 


sets  up  a  model  grid  over  a  designated  region  and  extracts  the  topographic 
and  bathymetnc  data  and  writes  a  TOPOG.DAT  file.  p^gi^pnic 


5.I.4.2 


Input  Conversion  Program 


The  input  conversion  program  extracts  the  various  required  fields  fi-om  the  NOGAPS  model  and 
outputs  all  necessary  variables  over  a  three-dimensional  region  at  each  time  in  direct  access 
ormat  The  output  files  are  wntten  for  each  time  covering  the  duration  of  the  required 
simulation.  The  naming  convention  for  the  files  is  NOGAPS.OOO  for  the  zero*  hour  fanalvciO 
data  file,  NOGAPS.006  for  the  6*  hour  data  file  and  so  on.  ^  ^  ^ 


5.L4.3 


The  Field  Generator 


me  field  generator  reads  the  TOPOG.DAT  file  and  inputs  the  atmospheric  data  files  sequentially 
(either  derived  from  the  NOOAPS  model  or  a  coarse  resolution  simulation  of  MAPS)  arid 
honzontally  interpolates  the  fields  to  the  region  defined  by  the  TOPOG.DAT  parameters  The 
output  files  are  wntten  in  direct  access  format  to  AMSUBGRD.nnn,  where  “nnn”  is  a  three-diait 
integer  denoting  the  number  of  hours  from  initial  time  ° 

5. 1. 4. 4  The  Preprocessor 


A  T  MAPS.dat  parameter  file  to  obtain  the  vertical  levels  required  for 

the  MAPS  simulation  and  interpolates  the  AMSUBGRD  files  to  the  appropriate  vertical  levels 
here  are  two  versions  of  the  preprocessor,  AMPREPPL  and  AMPREPSL.  The  first  of  these 
assumes  input  on  pressure  levels  and  interpolates  the  AMSUBGRD.nnn  files  fi-om  pressure 
evels  to  sigma  levels,  the  second  of  these  assumes  input  on  sigma  levels  if  the  input  atmospheric 

if »  m  The  preprocessor  oulpuls  a 

senes  of  AMPREP.nnn  files  at  the  appropnate  time  interval  for  the  duration  of  the  run. 


5. 1.4.5 


Display  Modules 


There  are  options  for  displaying  MAPS  model  output.  The  first  of  these  is 
the  simulated  vanablcs  at  designated  levels  and  times  during  the  model 
second  TSPLOT  that  displays  the  time  varying  station  output. 


AMDISP  that  plots 
simulation  and  the 
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5.2  CSC  GCOM2D 

CSC  GCOM2D  is  a  depth-integrated,  barotropic  hydrodynamie  model  used  for  modeling  sea 
surfaee  height  and  mean  current  structure  in  coastal  regions  due  to  tides  and  atmospheric  forcing. 
It  is  economic  to  run  in  terms  of  computational  overhead  and  data  initialization  requirements. 

5.2.1  Constraints  and  Limitations 


GCOM2D/3D  has  been  designed  such  that  total  grid  arrays  in  the  horizontal  direction  are  no  larger 
than  40000  points  (for  example,  200x200). 

5.2.2  Logic  and  Basic  Equations 

The  equations  solved  are  the  shallow  water  equations  and  are  presented  below; 


m  dP 


dx  a.v 


du  „at/ 


-m  U  —  +  V— 
dx  dy 


^sx  ^bx 


dx  dy 


dv  dC  m  dP 


p(o  V 

"’■  a.v  av 


(,,dV  ^dV 

-m  U - +V - 

.  a.v  dy 


a  V  d  fVH 


dx[  m  j  ai 


(5.2.1) 


(5.2.2) 


(5.2.3) 


where  U  and  V  arc  the  deptli  averaged  currents  in  the  x  and  y  directions  respectively,  H  is  the  total 
depth,  ^  is  the  surface  elevation,  /  is  the  Coriolis  parameter,  g  is  the  acceleration  due  to  gravity,  P 

is  the  aunospheric  surface  pressure,  p^.  is  the  water  density,  u  is  the  coefficient  of  lateral  eddy 
diffusion  and  has  a  value  of  0.2,  and  r^,,,  the  bottom  frictional  stress  in  tlie  x  and  v  directions, 
respectively.  Sxx,  S^y,  ^^  and  5,.,  represent  tlie  wave  radiation  stresses. 

The  equations  have  been  fomiulated  on  a  Cartesian  grid  due  to  the  relative  case  of  coding. 
However,  the  fonnulation  is  generalized  to  incorporate  different  map  projections  through  the 
appropriate  specification  of  the  map  factor,  m,  a  scaling  factor  dependent  on  the  chosen  map 
projection  of  the  model  grid.  In  accordance  with  this  generality,  the  Coriolis  parameter  varies  with 
latitude. 

The  surface  wind  .stress  components  arc  computed  using  the  quadratic  relationship; 
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=  CDPaK\u^,  =  C^pjujv,, 

(5.2.4) 

where  |u  J  =  (u]  +v^'f"  Ua  and  are  the  horizontal  components  of  \\^ind  velocity  at  anemometer 

height,  IS  the  density  of  air,  and  Cd  is  the  drag  coefficient  based  on  Smith  and  Banke  (1 975)  and 
expressed  as  follows; 

Q  =[0.63  +  0.066|uj]xl0-^’  |uj<25ms-'; 

=[2.28  +  0.033(|uJ-25)]xl0-^-  |uj>25ms-’. 


^e  bottom  stress  is  represented  by  a  Manning’s  n  depth-dependent  fiiction  relation 
Signell  and  Butman  (1992): 


following 


(5.2.6) 


where  n  has  the  value  0.03.  This  formulation  ensures  that  the  drag  coefficient  increases  with 
decreasing  water  depth  and  is  applied  to  water  depths  greater  than  1  m.  In  extremely  shallow  water 
and  over  land  points  that  become  inundated,  drag  coefficients  can  be  specified  at  each  grid-point 
according  to  the  terrain  type.  ^  ^ 

5.2.2. 1  Solution  Procedure 

Equations  (5.2  1)  to  (5.2.3)  are  solved  using  a  split-explicit  finite-difference  scheme  on  an 
rakawa-C  ^d  (Messinger  and  Arakawa,  1976)  as  shown  in  Figure  5.2.2- 1  and  described  in 
Hubbert  et  al.  (1990).  The  continuity  equation  and  the  gravity  wave  and  Coriolis  terms  in  the 
momentum  equations  are  solved  on  the  shortest  time  step,  known  as  the  adjustment  step,  using 
the  forward-backward  method.  The  non-linear  advective  terms  are  solved  on  an  intermediate 
advective  time  step  using  the  two-time-level  method  of  Miller  and  Pearce  (1974).  Finally  on  the 
longest  time  step,  the  so-called  physics  step,  the  surface  wind  stress,  bottom  friction  stress  and 
atmospheric  pressure  terms  are  solved  using  a  backward-implicit  method.  This  approach  is 
extremely  efficient  in  oceanographic  models  with  free  surfaces  because  of  the  large  disparity 
between  adveelive  speeds  and  gravity-wave  phase  speeds  in  deep  w^ater. 
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Figure  5.2. 2-1:  The  horizontal  grid  structure  used  in  GCOM2D. 

(a)  Adjustment  step: 

By  defining 

m'.  =(m^.  +n7^.^,)/2 
//:=(//,+//,,„)/2 
H:  =(//,+//, ,,J/2 


and  noting  that  the  grid  spacing  at  the  standard  latitude  of  the  map  projection  is  represented  by 
Av  =  Av  =  Ar ,  the  adjustment  step  is  accomplished  by  the  following; 
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topKdSSpr'' 


y  y 


u*.=uuJjv:-m';JK 


(5.2.8) 


y  y  r  "  dy 


(5.2.9) 


These  equations  may  be  rewritten  directly  in  the  explicit  form  to  obtain  the  intermediate 
solutions  U  and  V  ; 


*  1  ^ 
U..  = - ; - 

y  (1  +  /'A/')  ij 


(5.2.10) 


*  1  f  r)/'  I 


Where 


(5.2.11) 


Equati^s  (5.2.^,  (5.2.10)  and  (5.2.1 1)  arc  solved  times  so  that  the  adjustment  step  advances  the 
intermediate  solution  to  the  values  U*  and  V*  at  time  t  +  A/„ ,  where  Ar^  =  N  At . 

(b)  Advective  step: 

pLme'TiTTS'  advective  step  is  the  two-time-level  method  of  Miller  and 

1  earce  (1974).  This  scheme  alternates  the  Euler  and  Euler-backward  schemes  at  odd  and  even 

un^frulrro'*"^',  *i!  ""  amplification  factor  that  is  almost  exactlv 

unit>  for  the  Courant  numbers  that  are  found  in  ocean  models  (see  Hubbert  et  al.  1991). 

On  the  odd  advective  time-steps,  the  updated  intennediatc  solutions  6'**  and  1'**  are 


U..  =(u;+At^F(U-l 


I'..  =  (V'  +AtG(V') 

IJ  y  a  \  if  /  a 


(5.2.12) 


(5.2.13) 


40 


PSI  Technical  Report  SSC-005-00 


PCTides  SDD 


where  the  advective  operators  F  and  G  are  defined  by 

+ v;  w:,»  +£/■)+ (v;.,  +  c,-,  w,  + 1/,;., » 


*  ni- 

Gi'-y- ) = vju  +  -  n- ) + + uUj.,  xv;  -  k„  ) + 


On  the  e\en  advective  time-steps,  the  two-step  Euler-backward  scheme  is  applied  instead  of 
equations  (5.2.12)  and  (5.2.13).  On  the  second  iteration,  the  values  of  U  and  V  obtained  after  the 
first  iteration  are  substituted  into  the  advective  operators  F  and  G  to  obtain  the  final  values  of  U** 
and  V** 


(c)  Physics  step: 

The  adjustment  and  advective  integration  cycle  is  carried  out  Np  times  so  that  the  interim 
solution  is  U**  and  V**  at  time  t  +  A/^,  where  =X^N^At.  At  this  time,  the 

solution  cycle  is  completed  with  the  inclusion  of  the  physics  terms  using  a  numerical  technique 
similar  to  that  described  for  the  adjustment  step.  The  final  value  of  U  is  therefore 

* .  n  M 

1  ffhi 

(5.2.14) 


At 


•J  p. 


At . 


P. 


h;  h; 


J^u  IJ  ^  IJ  tj  '  ^  V^/  +  ly  ^  tj  ' 


Rewriting  in  explicit  form  gives 


t  +  A/  ..  A/  „ 
V ..  ^  =f/.. 

'  P. 


T  f^hi 

_ p  ^  \ 

//“  As  ^  ^ 


u  ^  U  ' 


Similarlv 


'>■  "  P, 


- 1 

_ 1 

/ 

1 

< 

--- — --(P-,-p;;) 

/ 

1+  "  ATf7  -+F 

( — 

1 _ 

/ 

//; 

The  numerical  scheme  therefore  is  split-explicit  and  consists  of  three  distinct  time  steps  At ,  A/„ 


and  A/ „  w  here  At  <  At  <  At. 

P  'I 
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5.2*2.2  Boundary  and  Initial  Conditions 


Boundaty  ^nditions  can  be  applied  in  a  range  of  ways  depending  on  the  type  of  process  beinn 
modeled.  Meteorological  forcing  is  applied  via  the  wind  stress  and  surface  p^sure  Sent 
m  equations  (5.2.1)  to  (5.2.3)  at  all  submerged  model  grid-points  in  the  compLtional  dSTn 

Tidal  and  meteorological  forcing  at  lateral  boundaries  is  achieved  by  specifying  the  incremental 
disp  acement  ^  the  water  surface  due  to  changes  in  tidal  height  and  atm^Xrfc  TSemi 
displacement.  The  lateral  boundaiy  conditions  are  applied  using  a  “one-way  nesting”  technique  The 
bound^  conditions  are  applied  to  the  appropriate  model  variable  <p  (representing  heil  or 
velocity)  with  decreasing  intensity  from  the  boundary  to  some  specifi^  number  of  Ldel  ^d 
points  (typically  around  1 2)  into  the  domain  according  to  the  following  equation;  ^ 


where  ^  is  the  model  predicted  value  and 
according  to  a  cosine  function  such  that 


(5.2.15) 

is  the  prescribed  boundary  value  and  oc  is  varied 


a=0.5(cos;r(l-«//j„^)+l),  n=\,n^^ 


(5.2.16) 


At  coasts  boundaries  and  along  riverbanks,  the  wetting  and  drying  of  grid  cells  is  accomnlished  via 
he  inundauon  a  gonthm  described  in  the  next  section.  On  ouSow,  a  rTdiation 
as  descnbed  in  Miller  and  Thoipe  (1981)  is  applied  to  the  velocity  field  to  prevent  th^b^d 
wave^nergy  within  the  numencal  domain,  while  on  inflow  boundaries,  a  zem-gradient  conlonl 


Under  any  prescnbed  boundaiy  forcing,  GCOM2D  is  initialized  by  setting  velocities  to  zero  and 
in  erpolates  global  Finite  Element  Solutions  versions  95. 1/2. 1  (FES  95. 1/2. 1)  (Shum  et  al  1 997) 
ekva  ion  field  to  the  model  grid.  It  is  customary  to  allocate  an  initial  spin-up  plrioTto  atw  die 

flin  T^e  snru^im '""n  throughout  the  computational 

without  windsand24  hrs  wl\Ss  "  hours 

5.2.2.3  Inundation  Algorithm 

The  ad^sment  of-ccasllmes  lo  account  for  flooding  and  draining  rakes  place  after  equations 
(5.2.  H5_ -3)  have  been  solved  as  desenbed  above.  The  coastal  boundary  is  confitnired  to  pass 
through  he  velocity  gnd-points  on  the  staggered  grid  in  a  stepwise  manner  such  tat  i,  Lssd 

hrough  the  (/points  in  the., .-direction  and  through  F points  in  then-direction.  The  velocities  on 
(he  boundaries  arc  assumed  lo  be  zero.  ciociiics  on 
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The  first  step  is  to  calculate  the  distance  in  the  x-  and  j^-directions  that  fluid  could  travel  in  a  time 
step  at  each  ^  grid-point  that  is  adjacent  to  the  coast.  The  depth-averaged  current  velocity  used 
in  this  calculation  is  taken  at  the  first  grid-point  on  the  seaward  side  of  the  ^  grid-point.  The 
travel  distance  is: 


AT”.  =  i 

l,j  ‘ 


Ax:j+Atu;'_,j,u>0 
AXl]'  +AtU"j,U  <  0 
AY.y'  +A/T;.”._,,r>o 
[Ar.j'+AfF”j,F<0 


(5.2.17) 

(5.2.18) 


where  AX"/  and  A?;."  '  are  the  distances  in  the  x  and  y  directions  that  the  water  traveled  in  the 

previous  time  step.  By  factoring  in  the  travel  time  of  the  fluid,  inland  grid  cells  are  prevented 
from  automatically  becoming  inundated  at  the  first  instant  the  water  level  at  the  coast  exceeds 
the  height  of  the  adjacent  dry  points.  In  equations  (5.2.17)  and  (5.2.18),  the  first  option  applies 
to  flooding  in  the  positive  x  or  y  direction  (or  draining  in  the  negative  .r  and  y  direction)  while  the 
second  option  applies  to  flooding  in  the  negative  x  and  y  direction  (or  draining  in  the  positive  x 
and  y  direction). 

The  testing  for  coastline  movement  proceeds  in  thex-  and  y-directions  separately.  If  the  height  of 
the  water  at  the  first  z-point  seaward  of  the  coastline  exceeds  the  topographic  height  at  the  first  z- 
point  landward  of  it,  and  the  accumulated  distance  traveled  in  the  gi\’en  direction  exceeds  the 
grid  separation,  then  a  new  sea-point  is  added  to  the  computational  domain.  The  velocity  at  the 
newly  acquired  velocity  point  is  extrapolated  from  adjacent  points  and  the  continuity  equation  is 
solved  to  obtain  the  water  depth  at  the  new  height  point.  Finally,  and  are  set  to  zero. 

The  procedure  for  draining  is  similar.  If  the  height  of  the  fluid  at  the  height  point  adjacent  the 
boundary'  drops  below  some  arbitrary  positive  height,  e  and  the  accumulated  distance  traveled 
by  the  fluid,  exceeds  the  grid  length,  then  draining  is  assumed  to  have  occurred.  The  height  grid- 
point  is  reclassified  as  dry  (i.e.  ^  =  0)  and  the  boundary  relocated  to  the  adjacent  wet  velocity 

grid-point.  Examples  of  the  performance  of  the  inundation  algorithm  can  be  found  in  Hubbert 
and  Mclnnes  ( 1 999a,  b) 

5. 2. 2. 4  Tidal  Data  Assimilation 

In  order  to  improve  the  simulation  of  tidal  forced  dynamics  a  facility  is  included  to  “nudge”  the 
model  solution  witli  tidal  height  predictions  at  tidal  stations  within  the  model  domain.  Global  tidal 
station  constituent  data  is  stored  in  ASCII  fomiat  in  the  file  TCANALS.DAT.  New  stations  can 
easily  be  added  by  adding  to  this  file  in  the  appropriate  fomiat. 

During  a  model  run  tlie  TCANALS.DAT  file  is  scanned  by  subroutine  TIDEOBS  during  its  first 
call  and  the  tidal  constituent  parameters  for  stations  witliin  the  model  domain  are  stored.  At  each 
model  time  step  the  tidal  height  is  predicted  at  each  station  and  used  to  “nudge”  the  model  solution. 
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The  nudging  method  is  based  on  denvmg  a  new  solution  at  grid  points  near  each  tidal  station  from  a 
wei^ted  combination  of  the  model  solution  and  the  station  sea  level  prediction.  The  weighting 
fonction  IS  calculated  from  the  product  of  the  parameter  WGTFAC  (default  value  =  0  50)  aS  th! 
half  cosine  function  (range  -  0  to  1)  used  for  nesting  (equation  5.2. 1 6).  The  weighting  of  the  station 
sea  level  prediction  goes  to  zero  at  a  defined  distance  from  each  station  determined  bv 
p^ameter  ZRADIUS  (default  value  =  40  kms).  The  values  of  WGTFAC  and  ZRADIUS  can  be 
changed  m  the  parameter  file  ASSIM.DAT  that  resides  in  the  WORK  directory.  Changes  to  these 
parameters  are  not  advised  however  unless  there  is  a  specific  reason. 


5.2.3 


Data  Elements 


5.2.3. 1  Commonly  Used  Parameters 
File:  2dparams.for 


Variable 

Type 

Description 

nk 

Integer 

Maximum  vertical  dimension. 

ni 

Integer 

Maximum  horizontal  dimension  (x-direction) 

M 

Integer 

Maximum  horizontal  dimension  tv-directinnt 

nipl 

Integer 

NI+1.  - - 

-njpl 

Integer 

"nj+i.  “  - - - - — 

5. 2. 3. 2  Subroutine  AD  VECT 

the  MAfrI  proCT^^^  advective  terms  of  the  momentum  and  continuity  equations.  It  is  called  by 


Calling  Sequence:  subroutine  advect(nx,  ny,dt) 

Data  Declaration:  real  dt 

integer  nx,  ny 


Arguments:  nx  Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 
Time  step. 

Common  Blocks:  2DCOMMON 
5. 2. 3. 3  Subroutine  ASTROZ 


This  routine  calculates  the  frequency  (in  cycles/liour),  tlic  astronomical 
the  amplitude  for  the  given  tidal  constituents.  It  is  called  by  TIDEOBS. 


arguments  (in  cycles),  and 
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Calling  Sequence:  subroutine  astroz  (Istc,  con,  freq,  indx,  utc,  \^c,  ftc,  ntemax,  neons,  nobs, 
hour,  xlat) 

Data  Declaration;  real  con,  freq,  utc,  vtc,  ftc,  hour,  xlat 

integer  Istc,  indx,  ntemax,  neons,  nobs 

Logical  unit  which  reads  the  standard  tidal 
constituent  data  from  STCONTHP.DAT. 

Tidal  constituents  array;  con  (ntemax,  nobs). 

Frequencies  of  tidal  constituents  calculated  in 
ASTROZ. 

Stored  look-up  table  indices  from  frequency  calculation. 
Nodal  correction  phase  calculated  in  ASTROZ. 

Astronomical  argument  in  cycles  calculated  in  ASTROZ. 
Amplitude  of  tidal  constituents  calculated  in  ASTROZ. 
Maximum  number  of  tidal  constituents. 

Number  of  tidal  constituents  a\  ailable  at  each  tidal  station. 
Number  of  station  observations  available. 

Double  precision  hour  relative  to  time  origin  (01-01-1900). 
Station  latitude. 

Common  Blocks:  2DCOMMON 
5.2.3.4  Subroutine  ATMOS 


Arguments;  Istc 

con 

freq 

indx 

utc 

vtc 

ftc 

ntemax 

neons 

nobs 

hour 

xlat 


ATMOS  reads  the  meteorological  input  file  and  calculates  the  wind  stress  and  atmospheric  pressure. 
It  is  called  by  the  MAIN  program. 

Calling  Sequence:  subroutine  atmos(lu,  wdfile,  istart,  wdmax,  rhoa,  punits,  wtzone, 
ihour,  iminut,  iday,  imonth,  iyear,  nx,  ny,dt) 

Data  Declaration;  real  wdfile,  wdmax,  rhoa,  punits,  wtzone,  dt 
integer  nx,  ny, 

integer  istart,  ihour,  iminut,  iday,  imonth,  iyear 
logical  lu 


Arguments:  lu 

wdfile 

ist<ut 

wdmax 

rhoa 

punits 

wtzone 

iliour 

iminut 


Logical  unit  number  of  the  model  run  log  file. 

Name  of  meteorological  input  file. 

No.  of  hours  of  model  spin-up  time. 

Wind  speed  cut-off  in  wind  stress  fonnulation. 
Density  of  air. 

Multiplier  to  convert  pressure  from  pascals  into  mks 
units. 

Wind  time  zone,  in  real  hours  relative  to  GM'F. 

2  digit  integer  for  current  hour. 

2  digit  integer  for  current  minute. 
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iday 

imonth 

iyear 

nx 

ny 

dt 

Common  Blocks:  2DCOMMON 

5.2.3.5  Subroutine  CYCLE 


2  digit  integer  for  current  day. 

2  digit  integer  for  current  month. 
2  digit  integer  for  current  year. 
Array  dimensions  in  x-direction. 
Array  dimensions  in  y-direction. 
Time  step. 


airays  at  the  end  of  each  time  step.  It  is  called  by  MAIN 

Calling  Sequence:  subroutine  cycle(nx,  ny) 

Data  Declaration;  integer  nx,  ny 

Arguments:  nx  Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Common  Blocks:  2DCOMMON 

5. 2. 3. 6  Subroutine  ENERGY 


ENERGY  integrates  the  total  energy  of  the  system 
by  the  MAIN  program. 


as  a  diagnostic  performance  indicator.  It  is  called 


Calling  Sequence;  subroutine  energy(nx,  ny,egp) 

Data  Declaration,  real  egp 

integer  nx. 

Arguments;  nx 
ny 

egp 

Common  Blocks:  2DCOMMON 

5.2.3. 7  Subroutine  FIL TZ 


Array  dimensions  in  x-direction. 
Array  dimensions  in  y-direction. 
The  total  energy  of  the  system. 


FILTZ  applies  a  five-point  filter  to  the  sea  level 
by  the  .M.AIN  program. 


array  to  remove  two-grid-length  noise 


It  is  called 


Calling  Sequence:  subroutine  filtz(nx,  ny,  fltfac,  spv) 
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Data  Declaration:  real  fltfac,  spv 

integer  nx,  ny 


Arguments:  nx 
ny 

fltfac 

spv 


Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Inactive  argument. 

Special  value.  No  filtering  is  applied  when  height  = 
spv. 


Common  Blocks:  2DCOMMON 
Algorithm: 

See  section  5.7-  Filtering  Routines 
5.23.8  Subroutine  FILTUV 

FILTUV  applies  a  five-point  filter  to  the  velocity  arrays  to  remove  two-grid-length  noise.  It  is 
called  by  the  MAIN  program. 

Calling  Sequence:  subroutine  filtuv(nx,  ny,  fltfac,  spv) 

Data  Declaration:  real  fltfac,  spv 

integer  nx,  ny 


Arguments:  nx 
ny 

fltfac 

spv 


Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Inactive  argument. 

Special  value.  No  filtering  is  applied  when  height  = 
spy. 


Common  Blocks:  2DCOMMON 
5. 2. 3. 9  Subroutine  GRDSET 

GRDSET  sets  up  the  staggered  bathymetry,  coastlines  and  vertical  le\els  infonnation  (GCOM3D 
only).  It  is  called  by  tlie  MAIN  program. 

Calling  Sequence:  subroutine  grdsctfnx,  ny,dxyinin,  dxymax,  platO,  plongO,  platl,  plongl) 

Data  Declaration:  real  dxymin,  dxymax,  platO,  plongO,  platl,  plongl 
integer  nx,  ny 
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Arguments:  nx 
ny 

dxymin 

dxymax 

platO 

plongO 

platl 

piongl 


Array  dimensions  in  x-direction. 
Array  dimensions  in  y-direction. 


Value  of  minimum  grid  spacing  on  grid  for  lo^ 
Val  ue  of  maximum  grid  spacing  on  grid  for  loj 
Latitude  of  minimum  grid  spacing  (reporting 
variable). 

Longitude  of  minimum  grid  spacing  (reporting 
variable). 


Latitude  of  maximum  grid  spacing  (reporting 
variable). 

Longitude  of  maximum  grid  spacing  (reporting 
variable). 


Common  Blocks:  2DCOMMON 


5.2.3.10  S ubroutine  MASSBAL 


The  subroutine  monitors  the  mass  balance  throughout  the  model  integration  as 
performance  indicator. 

Calling  Sequence;  subroutine  massbal(nx,  ny,dt,  cum_sqkm,  iplot) 

Data  Declaration:  real  dt,  cum_sqkm 

integer  nx,  ny, iplot 

Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Model  time  step. 

Redundant  variable. 

Real-time  output  flag  (0=text  or  1  =  graphics  to 
screen). 


The  total  \  olume  of  the  computational  domain  at  any  time  is  given  by; 


Arguments:  nx 
ny 
dt 

cum__sqkm 

iplot 


Common  Blocks:  2DCOMMON 


nv  nx 


j-\  .=1 


(5.2.3.10a) 


The  mass  tlu.x  through  the  lateral  boundary  is  calculated  from  equation  (5.2.3. 
time  inteiv  al  dt.  The  total  mass  added  to  the  system  is  monitored  bv  comparin 
(5.2.3.10a)  and  (5.2.3.10b).  ' 


file. 

file. 


a  diagnostic 


10a)  over  a 
g  equations 
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Ki  =  'L^ojXdyx{H  +  ()„jXdl 

1=0 

ny 

+  ^dyx{H  +  .  xdt 

j=o 

nx 

+  'Zy,oXdyx{H  +  C\oXdt 

i=o 

tvc 

+  'L^,„xdyx(H  +  C),„xdt 


(5.2.3.10b) 


The  change  in  volume  ^/Fbetween  any  two  time  steps  n  and  n+1  is  gi%'en  by 


dVo,  =  KT' 


■  V" 

^ol 


(5.2.3.10c) 


5.2.3.11  Subroutine  PARAMS 

PARAMS  reads  the  parameter  file  GCOM.DAT  and  the  time  series  output  stations  file 
STATIONS.DAT.  It  then  determines  which  tidal  files  are  present  (if  required)  and  reads  the  initial 
wind  data  (if  required). 

Calling  Sequence:  subroutine  params(lu,  Iwd,  Hog,  bfile,  rsfile,  wdfile.  tefile,  wdmax, 

rhoa,  punits,  iwind,  itide,  iassim,  nstpts,  istart.  inund,  iplot,  idrift,  ibc, 
imass,  iflux,  nadv,  nphys,  ntemax,  ntes,  ntca.  nstmax,  sponge,  datum, 
fltfac,  dtfac,  wdt,  wtzone,  depmin,  depmax,  hcl.  hc2,  hc3,  hc4,  hc5,  ihour’ 
iminut,  iday,  imonth,  iyear,  hrout,  dtout,  nhours,  dhour,  dminut,'  dday! 
dmonth,  dyear,  dlat,  dlong,  sptime,  stname,  stnlat,  stnlng,  stdep,  thetap,  ip, 
jp,  stni,  stnj,  nstns,  region,  jproj,  proj,  height,  nx.  ny) 

Data  Declaration:  real  dt,  wdmax,  rhoa,  punits,  sponge,  datum,  fltfac. 

dtfac,  wdt,  wtzone,  depmin,  depmax,  hcl,  hc2,  hc3,  hc4,  hc5,  hrout,  dtout, 

dlat,  dlong,  sptime,  stnlat,  stnlng,  stdep,  thetap,  stni,  stnj, ,  proj,  height 

integer  nx,  ny, iwind,  itide,  iassim,  nstpts,  istart,  inund,  iplot,  idrift, 

ibc,  imass,iflux,  nadv,  nphys,  ntemax,  ntes,  ntca,  nstmax, 

ihour,  iminut,  iday,  imonth  iyear,  nhours,  ip.  jp,  nstns,  jproj,  dhour, 

dminut,  dday,  dmonth,  dyear, 

logical  lu,  Iwd,  Hog 

char  bfile,  rsfile,  wdfile,  tefile,  stname,  region 


Arguments:  lu 

Logical  unit  of  input  parameter  file. 

Iwd 

Logical  unit  of  wind  file  (ATMOS.DAT). 

Hog 

Logical  unit  of  output  log  file 

blilc 

Topography  file  name  (‘TOPOG.DAT'). 
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rsfile 

wdfile 

tefile 

wdmax 

rhoa 

punits 

iwind 

itide 

iassim 

nstpts 

istart 

inund 

iplot 

idrift 

ibc 

imass 

iflux 

nadv 

nphys 

ntemax 

ntes 

ntca 

nstmax 

sponge 

datum 

fltfac 

dtfac 

vvdt 

vvtzone 

depmin 

depmax 

hcl,  hc2,  hc3, 
hc4,  hc5 

iliour,  iminut,  iday 

imonth,  iyear 

lirout 

dtout 

nliours 


Restart  file  name  (‘GCOM.RST’). 

Atmospheric  input  file  name  (‘ATMOS.DAT’). 
Tidal  constituent  file  name  (‘M2.DAT,  S2.DAT’, 
etc.). 

Wind  speed  cut-off  in  wind  stress  formulation. 
Density  of  air. 

Multiplier  to  convert  pressure  fi-om  pascals  into  mks 
units. 

Wind  forcing  flag  (0=no  wind  forcing,  1=  wind 
forcing). 

Tide  forcing  flag  (0=not  tidal  forcing,  1=  tidal 
forcing,  2=  assimilate  data). 

Tidal  assimilation  flag  (1=  assimilated  tidal 
information). 

Number  of  points  for  nesting  in  (set  internally). 
Number  of  points  for  nesting  in  (set  internally). 
Inundation  flag  (l=call  inundation  algorithm). 

Screen  or  text  flag. 

Inactive  flag. 

Boundaiy  condition  flag  (set  to  5). 

Flag  for  calculating  mass  balance  (set  to  0). 

Inactive  flag. 

Number  of  gravity  wave  steps  before  advection  step. 
Number  of  gravity  wave  steps  before  physics  step. 
Maximum  number  of  tidal  constiments. 

Number  of  tidal  constituents  available  at  each  tidal 
Station. 

Inactive. 

Maximum  number  of  stations. 

Inactive. 

Inactive. 

Inactive. 

CFL  reduction  factor  (internally  set  0.95). 
time  stepping. 

Wind  time  zone,  in  real  hours  relative  to  GMT. 
Minimum  allowable  depth  in  model  (set  in 
PARAM). 

Maximum  allowable  depth  in  model  (set  in 
PARAM). 

Heights  tliat  define  color  bands  for  screen  graphics. 

Initialization  time  of  model. 

Initial  output  time  (set  at  0). 

Time-step  of  GCOM.OUT  file, 
total  number  of  hours  of  model  run. 

50 


PSI  Technical  Report  SSC-005-(X) 


PCTides  SDD 


dhour,dniinut, 

dday,  dmonth, 

dyear,  dlat,  dlong 

sptime 

stname 

stnlat 

stnlng 

stdep 

thetap 

ip 

jP 

stni 

stnj 

nstns 

region 

jproj 

proj 

height 

nx 

ny 


Inactive. 

Inactive. 

Inactive. 

Inactive. 

Inactive. 

Inactive. 

angle  of  the  grid  to  the  lat-long  grid  at  a  designated 
point. 

Integer  i-grid  points  value  of  output  stations. 
Integer  j-grid  points  value  of  output  stations. 

Real  i-grid  point  value  of  output  stations. 

Real  j-grid  point  value  of  output  stations. 

Number  of  stations 

80  character  title  from  the  TOPOG.DAT  file. 

Flag  defining  Lambert-conformal  projection. 

Array  of  Lambert-conformal  map  projections. 
Array  of  topography. 

Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 


Common  Blocks:  2DCOMMON 

5.2.3.12  Subroutine  OMOUT 

Calling  Sequence;  subroutine  omout(lu,  dfile,  region,  hour,  minut,  da\’,  month,  year,  tstep, 
wtzone,  ntes,  iwind,  ntemp,  punits,  nx,  ny,  nz) 

Data  Declaration:  real  tstep,  wtzone,  punits 

integer  nx,  ny,  nz,  hour,  minut,  day,  month  year,  ntes,  ntemp,  iwind 

logical  lu 

char  region,  dfile 


Arguments: 


lu 

Logical  unit  number  of  output  file. 

dfile 

Name  of  output  file  (GCOM.D.AT). 

region 

80  character  title  from  the  TOPOG.DAT  file. 

hour 

Integer  time  of  current  time  in  model  integration. 

minut,  day, 
month,  year 

Time  of  output. 

tstep 

Model  time  step. 

wt/one 

Wind  time  zone,  in  real  hours  relative  to  GMT. 

ntes 

Number  of  tidal  constituents  available  at  each  tidal 

iwind 

station. 

Wind  forcing  flag  (0=no  wind  forcing,  1  =  wind 
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forcing). 

ntemp 

Inactive. 

punits 

Multiplier  to  convert  pressure  fi-om  pascals  into  mks 
units. 

nx 

Array  dimensions  in  x-direction. 

ny 

Array  dimensions  in  y-direction. 

nz 

Array  dimensions  in  z-direction  (=1  for  GCOM2D). 

Common  Blocks:  2DCOMMON 
5.2.3. 13  Subroutine  PHYSIC 


calculates  the  adjustment  to  the  momentum  equations  due  to  surface  and  bottom  stress  and 
atmosphenc  pressure  gradient.  It  is  called  by  the  MAIN  program.  ^ 


Calling  Sequence;  subroutine  physic(nx,  ny,  iwind,  dt,  chev,  rhobar) 

Data  Declaration:  real  dt,  chev,  rhobar 

integer  nx,  ny,  iwind 


Arguments:  nx 

Array  dimensions  in  x-direction. 

ny 

Array  dimensions  in  y-direction. 

iwind 

Wind  forcing  flag  (0=no  wind  forcing,  1=  wind 
forcing). 

dt 

Time  step. 

chev 

1  f 

Coefficient  of  horizontal  eddy  \iscosity. 

rhobar 

Common  Blocks:  2DCOMMON 

Mean  density  of  water. 

5.2.3.14  Subroutine  TIDEOBS 


This  sub^routine  opens  the  tidal  station  data  file  TCANALS.DAT,  and  determines  which  stations  are 

Zdel  T  II  carries  out  tidal  height  predictions  at  each  station  and  nudges  the 

modeled  sea  levels.  It  is  called  by  the  MAIN  program.  ^ 

Calling  Sequence:  subroutine  tideobs(lobs,  Ithp,  tzone  Istc,  zem,  iassim,  tsdt,  hour  minut 
month,  year,  nx,  ny,dt)  ’ 


Data  Declaration;  real  dt,  tzone,  zerr,  tsdt 

integer  nx,  ny,  iassim,  hour,  minut,  month,  year 
logical  lobs,  Ithp,  Istc 


Arguments:  lobs 
Ithp 


Logical  unit  of  observation  file  input. 
Logical  unit  of  tidal  prediction  output. 
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Time  zone,  in  real  hours  relative  to  GMT. 
Logical  unit  of  L. 

Mean  error  between  model  and  tidal  stations. 
Assimilation  flag  (0=no  assimilation, 
l=assimilation). 

Time  series  output  time. 

Current  time  of  observations. 

Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Time  step. 


5.2.3.15  Subroutine  UVBAR 

UVBAR  calculates  w  at  v  grid  points  and  v  at  w  grid  points.  It  is  called  by  the  MAIN  program. 
Calling  Sequence:  subroutine  uvbar(nx,  ny) 

Data  Declaration:  integer  nx,  ny 

Arguments,  nx  Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 


tzone 

Istc 

zerr 

iassim 

tsdt 

hour,  minut 

day,  month,  year 

nx 

ny 

dt 

Common  Blocks:  2DCOMMON 


Common  Blocks:  2DCOMMON 
5.2.3.16  S ubroutin  e  UVBND 

This  subroutine  applies  boundary  conditions  to  velocity  components.  It  is  called  bv  the  MAIN 
program. 


Calling  Sequence:  subroutine  uvbnd(nx,  ny,  dt,  ibc,  nstpts) 

Data  Declaration:  real  dt 

integer  nx,  n\’.  ibc,  nstpts 


Arguments:  nx 

Array  dimensions  in  x-direction. 

ny 

Array  dimensions  in  y-direction. 

dt 

Model  time  step. 

ibc 

Boundary  condition  flag  (set  internally  to  5). 

nstpts 

Number  of  points  for  nesting  in  (set  internally) 

Common  Blocks:  2DCOMMON 
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5.2.5. 1 7  Subroutine  UVGRA  V 

UVGRAV  solves  the  gravity  wave  terms  in  the  equations  of  motion.  It  is  called  by  the  MAIN 
program. 

Calling  Sequence:  subroutine  uvgrav(nx,  ny,dt) 

Data  Declaration:  real  dt 

integer  nx,  ny 


Arguments,  nx  Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Model  time  step. 

Common  Blocks:  2DCOMMON 
5.2.3.18  Subroutine  ZBOUND 

ZBOUND  inputs  nesting  heights  on  boundaries  due  to  tides  (from  S2.DAT,  M2.DAT,  etc)  coarse 

mesh  model  data  (from  GCOM.NST)  or  inverse  barometer  effects  (from  ATMOS.DAT)  It  fr  called 
by  the  MAIN  program.  ’ 

Calling  Sequence:  subroutine  zbound(nx,  ny,iwind,  ipress,  itide,  nstpts.  ibc,  dt,  rhobar,  pbar) 

Data  Declaration:  real  dt,  rhobar,  pbar 

integer  nx,  ny,  iwind,  ipress,  itide,  nstpts,  ibc 

Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Wind  forcing  flag  (0=no  wind  forcing,  1=  wind 
forcing). 

Pressure  flag  (0=  no  atm.  pressure,  so  set  to  pbar) 

Tide  forcing  flag  (0=  no  tidal  forcing,  I  =  tidal 
forcing,  2=  data  assimilation). 

Number  of  points  for  nesting  in  (set  internally). 

Boundary  condition  flag. 

Model  time  step. 

Mean  density  of  water. 

Mean  atmospheric  pressure. 

Common  Blocks:  2DCOMMON 


Arguments:  nx 
ny 

iwind 

ipress 

itide 

nstpts 

ibc 

dt 

rhobar 

pbar 
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5.2.3.19  Subroutine  ZSOL  VE 

This  routine  solves  the  continuity  equation.  It  is  called  by  the  MAIN  program. 

Calling  Sequence:  subroutine  zsolve(nx,  ny,dt) 

Data  Declaration:  real  dt 

integer  nx,  ny 


Arguments:  nx  Array  dimensions  in  x-direction. 

ny  Array  dimensions  in  y-direction. 

dt  Model  time  step. 

Common  Blocks:  2DCOMMON 

5.2.3.20  Subroutine  ZVNEST 

ZVNEST  interpolates  the  coarse  grid  GCOM.OUT  to  fine  grid  GCOM.NST.  It  is  called  by  the 
MAIN  program. 

Calling  Sequence:  subroutine  zvnest(lu,  cgfile,nx,  ny,hour  minut,  day,  month,  year,  dt, 
spv,  ntesnst) 

Data  Declaration:  real  dt,spv 

integer  nx,  ny,  hour,  minut,  day,  month,  year,  ntesnst 
char  cgfile 
logical  lu 

Arguments:  lu 

cgfile 
nx 
ny 

hour,  minut,  day, 
montli,  year 
dt 
spv 

ntesnst 


Logical  unit  number  of  coarse  grid  file  interpolated  to 
fine  grid. 

Coarse  grid  file  name. 

Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Nesting  time. 

Model  time  step. 

Special  Value. 

Number  of  tidal  constituents  in  nesting  file. 


Common  Blocks:  2DCOMMON 
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5.2.3.21  Subroutine  COASTS 

prorar  ®  boundan-.  It  is  called  by  the  MAIN 

CaUing  Sequence:  subroutine  coasts(iix.  ny,  ibf,  dt,  depmin,  depmax,  iplot,  inund) 

Data  Declaration:  real  dt,  depmin,  depmax 

integer  nx,  ny,  ibf,  iplot,  inund 


Arguments:  nx 
ny 
ibf 
dt 

depmin 

depmax 

iplot 

inund 

Common  Blocks:  2DCOMMON 


Array  dimensions  in  x-direction. 

Array  dimensions  in  y-direction. 

Bathymetry  filtering  flag  (set  internally  to  ‘on’). 
Model  time  step. 

Minimum  allowable  depth  in  model  (set  in 
PARAM). 

Maximum  allowable  depth  in  model  (set  in 
PARAM). 

Real-time  output  flag  (0=text  or  1  =  graphics  to 
screen) 

Inundation  flag  (1=  call  inundation  algorithm). 
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5.2.4  The  G  COM  System 

In  addition  to  GCOM2D,  several  other  modules  and  programs  are  necessary  to  complete  the 
modeling  system.  These  include  programs  for  setting  up  model  grids  over  specified  geographic 
regions,  extraction  and  conversion  programs  to  provide  atmospheric  initial  and  boundary 
conditions  in  suitable  format  for  GCOM2D  and  display  programs  for  viewing  model  output. 
These  additional  programs  are  summarized  in  Figure  5.2.4-1  while  the  key  additional  elements 
are  described  in  the  following  sections. 


GCOM  running  environment 


Figure  5. 2.4-1:  Flowchart  illustrating  the  running  of  the  GCOM2D/3D  system. 
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5.2.4. 1  Field  Generator 

The  field  generator  reads  the  TOPOG.DAT  file,  inputs  the  atmospheric  data  files  (MAPS.OOO 
etc),  and  horizontally  interpolates  the  fields  to  the  region  defined  by  the  TOPOG  DAT 
parameters.  The  output  fields  are  written  in  sequential  file  format  to  ATMOS.DAT. 

5.2. 4. 2  Tide  Model 


The  tide  model  reads  the  global  30-minute  tidal  constituent  files  (GLBM2.30M  GLBS2  3niu 
GLBO1.30M,  GLBK1.30M  etc.)  derived  from  the  global  model.  The  Finite  Element  Solutions 
95.1.21  (FES  95.1/2^1)  (Shum  et  al.,  1997).  It  then  interpolates  the  data  for  each  constituent  to 
the  region  defined  by  the  TOPOG.DAT  parameters.  The  output  files  are  named  M2.DAT 

S2.DAT,  (Dl.DAT,  Kl.DAT  etc.  These  values  may  be  displayed  in  the  system  output  model  (see 
sect.  5. 2.4. 4). 


5.2.4.3 


Hurricane  Model 


The  primary  requirement  for  modeling  storm  surges  is  accurate  surface  wind  and  atmospheric 
pressure  fields.  There  are  usually  insufficient  data  to  allow  a  direct  analysis  of  the  central  region  of 
rnost  tropical  cyclones  and  currently  available  numerical  weather  prediction  models  do  not 
adequately  represent  the  small-scale  features.  Surface  wind  speeds  and  pressures  are  therefore 
denved  by  means  of  a  numerical  tropical  cyclone  model  based  on  the  empirical  relationships 
presented  by  Holland  (1980)  and  similar  to  the  model  described  in  Hubbert  et  al  (1991)  This 
method  of  deriving  tropical  cyclone  wind  and  pressure  fields  has  been  tested  extensively  and  has 
also  been  used  extensively  in  forecasting  offices  both  in  Australia  and  elsewhere  and  has  been  found 
to  provide  a  good  representation  of  the  wind  fields  in  the  vicinity  of  a  cyclone 


The  pressure  P  (hPa)  at  radius  r  is  derived  as  follows: 

P  =  Pc+iP„-Pc)e~^'-'^^' 


(5.2.19) 


where  is  tlie  central  pressure,  P„  is  the  environmental  pressure  (the  climatological  mean  for  the 
region  and  month),  /-„  is  the  radius  of  maximum  winds  and  b  provides  a  scaling  on  the  profile  shape 
The  parameter  6  IS  empirically  defined  by  ‘ 

6  =  1.5  +  (980-PJ/120  (5.2  20) 

The  symmetric,  gradient-level  azimuthal  wind  component  is  estimated  bv 


v=  r\f\ 

.  I  '•  j  P  "  4  J  2 

where  p  is  the  air  density  and  /"  is  tlie  Coriolis  parameter. 


(5.2.21) 


A  lirst  order  asymmctiy  is  included  by  adding  the  cyclone  translation  to  the  symmetric  field  and 
rotating  the  field  so  that  the  maximum  wind  is  70"  to  the  left  (right  in  the  Nonliem  Hemisphere)  of 
the  direction  ot  cyclone  motion.  Tlie  radial  wind  field  is  constructed  by  rotating  the  flow  to  a 
constant  inflow  angle  ol  25  outside  the  radius  of  maximum  w'inds. 
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It  is  important  to  note  that  the  tropical  cyclone  model  is  not  expected  to  represent  the  full  field  of 
synoptic  scale  features  with  a  high  degree  of  accuracy.  The  critical  aspect  for  storm  surge  and  wave 
forcing  is  that  the  model  parameterizes  the  mesoscale  forcing  in  the  vicinity  of  the  maximum  winds 
reasonably  well. 

5. 2. 4. 4  Display  Modules 

There  are  two  options  for  displaying  GCOM2D  output.  The  first  of  these  is  OMDISP  that  plots 
the  simulated  variables  at  designated  levels  and  times  during  the  model  simulation  and  the 
second,  TSPLOT  that  displays  the  time  varying  station  output  variables  (sea  level  and  currents). 
In  addition,  the  bathymetry  may  be  displayed  in  a  2-dimensional  form  using  BATHPLOT  or  in  a 
3-dimensional  form  using  PLOT3D.  Boundary  conditions  derived  fi'om  FES95. 1/2.1  (Shum  et 
al.,  1997)  can  be  displayed  with/without  IHO  tide  station  data  using  PLOTIDE. 

The  organization  of  subroutines  and  functions  within  the  GCOM2D/3D  model  is  illustrated  in 
Figure  5. 2. 4-2. 
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Fioure  5. 2.4-2:  continued. 
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GCOM2D/3D 
Main  Program 
(continued) 


Subroutines 


Subroutines/Functions 


JSPONGE><j: 


SETWGT 


Calculate  coefficient  of 
bottom  friction 


Adjust  bottom  friction 
coefficients 

I 

Output  initial  fields 


Write  initi  J  conditions 
to  log  file  , 


OPNOUT 


SEQOUT 


ISTART>0 


CYCLE 

Read  GCOM,DAT  and  set  flags 


Figure  5.2.4-2:  continued 
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GCOM2D/3D 
Main  Program 
(continued) 


1PLOT=0 


Write  label 
headings 


IPLOT=l 


Setup  2D 
plotting  to 
screen 


Subroutines 


get  setup  time 


BGNPLT 

Initialise  plotting  routines 


Figure  5.2.4-2:  continued. 


Subr  0  utines/F  unctions 
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ime-steppmg  Loop 
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GCOM2I>^I) 
Main  Program 
(continued) 

-►i 


Commence  timestepping 
Yes 


No 


No 


Update  boundary  z 


Yes 


Yes 


Subroutines 


Subroutines/Functions 


ADATE 

OMOUT 

ADATE 

Z] 

ATMOS 

Derive  time  dependent  wind 
stress  and  atmospheric  pressure 


CYCLE 

Cycle  time-dependent  arrays 


ZSOLVE 
Update  2 


thpred 


ZBOUND 


TIDEOBS 


Figure  5.2.4-2;  continued. 


ATIME 


GRDPRJ 


SEQINP 


STRESS 


Yes 

ZVNEST 

Read  residuals  along  open 
boundaries 

h- 

ATIME 

i 

■i 

SEQINP 

GRDPRJ 

ASTROZ 

ATIME 

U| 

DAT!  ME 
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GCOM2D/3D 
Main  Program 
(continued) 

Update  U  and  V 
(a)  Gravity  wave  step 


(b)  Advective  step 


^NADV>0>v 

INSTEP,NADV|=^ 


g'  (c)  Physics  step 


PSwSPINUPv 


vISTEP.NFILTp 


Subroutines 


WSOLVE  tCCOM^n  only 


EDVIS  fCCOM.tn  only 


FILTZ 

Apply  2-grid-Iength  filter  to  Z 
FILTUV 

Apply  2-grid-Icngth  filter  to  U,V 


Subroutines/Functions 


HnLT3 
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GCGM2D/3D 
Main  Program 
(continued) 


'HOURS>' 

vSPINUPv 


^STEP^OU 


Update  screen  graphics 


^STEPXTSOU 


- 4000  I 

4100"^ - 

x^Tiie  output  to  log  file 


Figure  5.2A-2:  continued. 


Subroutines 


Subroutines/Functions 
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5.3  CSC  GCOM3D 


GCOM3D  is  a  3-dimensional  z-coordinate  hydrodynamic  model  that  is  used  for  coastal  ocean 
modeling  where  current  structure  with  depth  is  required.  It  can  be  run  in  either  barotropic  mode 
(no  thermal  and  density  variation)  or  baroclinic  mode  (temperature  and  salinity  are  carried  as 
prognostic  model  variables).  However  only  the  barotropic  mode  is  supplied  with  this 
installation. 

5.3.1  Constraints  and  Limitations 


See  Section  5.2.1. 

5.3.2  Logic  and  Basic  Equations 


For  folly  three-dimensional  incompressible  flow  the  equations  of  momenrnm  conservation  are: 


dl 


I 


I  ^ 

— 


Pk 


\ 


ax  ay  az 


(5.3.1) 


dr 


k  '^k 


’'‘■^-yy<-’'kV-i<'^k''k> 


1 


dP. 


Pk 


pi 


dx 


d 

dy 


+—tP'  +- 


-T- 


(5.3.2) 


dll 


dr 


1 


dP. 


Pi.  Pi 


dr 


^k  + 


d  ,,  d  ,  ^ 


J 


(5.3.3) 


In  coastal  ocean  now.s,  there  is  little  vertical  motion  and  so  the  \ertical  acceleration  and  the 
vertical  gradients  of  the  stress  terms  arc  negligible  compared  to  g  so  equation  (5.3.3)  reduces  to 
the  hydrostatic  equation, 
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(5.3.4) 


The  equation  of  continuity  for  an  incompressible  fluid  is’ 


(5.3.5) 


In  order  to  close  this  set  of  equations,  the  stress  terms  must  be  expressed  in  terms  of  known 
vanables.  The  honzontal  shearing  stresses  are  usually  written  in  Boussinesq  form  as; 


T-=pK„ 


dU 
dx  ’ 


dy 


dz  ’ 
dV 
dx  ’ 

dy 


=pK 

r^=pKH 


(5.3.6) 


where  K  K{x,y,z,t)  is  the  coefficient  of  vertical  eddy  viscosity  and  is  the  coefficient  of 

honzontal  eddy  viscosity  and  is  assumed  constant.  Equations  (5.3.1)-(  5.3.6)  form  the  basis  for 
modeling  the  vertical  and  honzontal  variation  of  the  ocean  dynamics. 

iTJnn.T’  ^  ^  ^  the  fluid  and  there  are 

vanous  levels  of  approximation  which  can  be  adopted.  In  the  barotropic  version  of  GCOM3D  it 

IS  assumed  that  =  K(U,  V,  z)  =  C,  [{ff  y  {ff  J.  This  is  the  formulation  adopted  by  Leendertsee 

(1970).  Another  alternative  is  to  assume  that  K  is  a  constant  however,  this  approach  is  not  found 

0  ^ye  realistic  vertical  mixing.  A  more  sophisticated  approach  is  to  solve  for  K  using  the 

turbulent  closure  models  described  in  Mellor  and  Yamada  (1974,  1982)  The  snecification  of 

surface  stress  and  bo«om  stress  in  GCOM3D  is  identical  that ’of  GCOM2D  (s"  “ZLns 


5.3.3  Data  Elements 


5.3.3. 1  Commonly  Used  Parameters 


File:  3dparams.for 


Variable 

Type 

Description 

nk 

Integer 

Maximum  vertical  dimension 

ni 

Integer 

Maximum  horizontal  dimension  fx-direction) 

nj 

Integer 

Maximum  honzontal  dimension  (y-direction) 
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Subroutines 

5.3.3.2  Subroutine  EDVIS 

EDVIS  calculates  the  coefficient  of  vertical  eddy  viscosity  at  each  grid  point. 

Calling  Sequence;  subroutine  edvis  (iev,  iwind,  itide,  ny,  zk,  cevl,  cev2,  wk,  s,  eg,  cbf) 

Data  Declaration:  real  zk,  cevl,  cev2,  wk,  s,  eg,  cbf 
integer  iev,  iwind,  itide,  ny 

Eddy  viscosity  formulation  flag  (internally  set). 

Wind  forcing  flag(0=  no  wind  forcing,  1=  wind 
forcing). 

Tide  forcing  flag  (0=  no  tidal  forcing,  1=  tidal  forcing,  2= 
data  assimilation). 

Number  of  grid  points  in  the  y-direction. 

Constant  eddy  viscosity  (not  used). 

Inactive. 

Coefficient  of  vertical  eddy  viscosity. 

Inactive. 

Inactive. 

Inactive. 

Inactive. 

Common  Blocks:  3DCOMMON 
5.3.3.3  Subroutine  WSOLVE 

WSOLVE  calculates  the  vertical  current  velocity. 

Calling  Sequence;  subroutine  wsol\'e(ny) 

Data  Declaration;  integer  ny 

Arguments:  ny  Number  of  grid  points  in  y-direction. 

Common  Blocks;  3DCOMMON 

Library  Programs  and  Data  Files 


Arguments:  iev 
iwind 

itide 

ny 

zk 

cevl 

cev2 

wk 

s 

eg 

cbf 


Libraries  arc  used  to  store  frequently  u.scd  subroutines  that  arc  common  to  various  PCTides 
programs.  This  section  documents  the  commonly  u.sed  sets  of  library  routines. 
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5.4  Input/Output  Routines 


(a)  Direct  Access  Files 

The  I^S  ptem  uses  a  direct  access  file  format  where  all  model  variables  at  all  levels  are 
stored  at  each  given  time.  There  are  various  subroutines  associated  with  the  reading  and  writing 

of  direct  access  files  and  these  are  stored  in  the  library  AMSUBS.FOR.  A  brief  description  of 
these  routines  follows.  ui 


5.4.1  Subroutine  DAOPEN 


This  routine  opens  the  direct  access  file  connected  to  unit  7m’. 

Calling  Sequence:  subroutine  daopen  (luin,fiiame,nxy) 

Data  Declaration:  logical  luin 

integer  nxy 
char  fiiame 

Arguments:  luin  Logical  unit  number  of  input  file  name. 

fhame  Character  variable  containing  file  name, 

nxy  Total  record  length. 


5. 4. 2  Subroutine  DA  CLOS 

This  subroutine  closes  the  direct  access  file  connected  to  unit  ‘LU’. 


Calling  Sequence:  subroutine  dados  (lu) 

Data  Declaration:  logical  lu 

Arguments:  lu  Logical  unit  number  of  file  to  be  closed. 


5.4.3  Subroutine  DAPRMS 

working  elements  of  the  file  directory,  specifically  the  grid  projection 

Calling  Sequence:  subroutine  daprms(idir,rdir,cdir,isize,nx,ny,nzjproj.proj) 

Data  Declaration:  integer  idir,  isize,  nx,  ny,  nz,  jproj 
real  proj,  rdir 
char  edir 
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Arguments:  idir 

Integer  file  header  information  array. 

rdir 

Real  file  header  information  array. 

cdir 

Character  file  header  information  array. 

isize 

Record  length. 

nx 

The  dimension  in  the  x-direction. 

ny 

The  dimension  in  the  y-direction. 

nz 

The  dimension  in  the  vertical  direction. 

jproj 

Projection  parameter  flag. 

proj 

Element  array  of  projection  parameters. 

5.4.4  Subroutine  DADIR 

This  subroutine  sets  up  the  working  elements  of  the  file  directory. 


Calling  Sequence:  subroutine  dadir(idir,rdir,cdir,isize,slf,ns,mlf,nm,rlev,nl) 

Data  Declaration:  integer  idir,  isize,  ns,  mlf,  run,  nl 
real  rdir,  slf,  rlev 
char  cdir 


Arguments:  idir 

Integer  file  header  information  array. 

rdir 

Real  file  header  information  array. 

cdir 

Character  file  header  information  array. 

isize 

Record  length. 

slf 

Number  of  single  level  fields. 

ns 

Array  of  single  level  field  names. 

mlf 

Number  of  multi-level  fields. 

nm 

Array  of  multi-level  field  names. 

rlev 

Number  of  levels. 

nl 

Array  of  sigma  or  pressure  levels. 

5.4.5  Subroutine  DALEVS 

This  routine  extracts  the  levels  information  from  the  directory. 

Calling  Sequence:  subroutine  dalcvs(idir,rdir,cdir,isizc,rlev,nl) 

Data  Declaration:  integer  idir,  isize,  nl 
real  rdir,  rlev 
char  cdir 

Arguments:  Sec  DADIR 
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5. 4. 6  Subroutine  DARDDR 

This  routine  obtains  the  directory  from  file  record  1. 

Calling  Sequence:  subroutine  darddr  (lu,idir,rdir,cdir,isize) 

Data  Declaration:  logical  lu 

integer  idir,  isize 
real  rdir 
char  cdir 

Arguments:  See  DADIR 

5. 4. 7  Subroutine  DA  WRDR 

The  purpose  of  this  subroutine  is  to  write  the  directory  to  file  record  1 . 

Calling  Sequence:  subroutine  dawrdr  (lu,idir,rdir,cdir, isize) 

Data  Declaration:  logical  lu 
real  rdir 

integer  idir,isize 
char  cdir 

Arguments:  See  DADIR 

5. 4. 8  Subroutine  DAINP 


This  routine  reads  the  field  as  specified  by  name  and  level.  If  IRET=0,  the  data  was  extracted 

successftilly,  IRET-1,  the  field  or  level  was  invalid  or  not  present.  IRET=2,  tliere  was  an  error 
while  reading  the  file. 


Calling  Sequence:  subroutine  dainp(lu,idir,rdir,cdir,name,rlev,array,ni,nx,ny,iret) 

Data  Declaration:  logical  lu 

real  rdir,  rlev,  array 
integer  idir,  ni,  nx,  ny,  iret 
char  cdir,  name 
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Arguments:  lu 
idir 
rdir 
cdir 
name 
riev 
array 
ni 
nx 
ny 
iret 


Logical  unit  number  of  input  file. 

Integer  file  header  information  array. 

Real  file  header  information  array. 

Character  file  header  information  array. 

4  Character  string  containing  name  of  field  to  be  read  in. 

1- D  array  containing  vertical  levels. 

2- D  array  to  be  read  in. 

The  first  dimension  of  ARRAY. 

The  number  of  elements  in  the  1**  dimension  of  ARRAY. 
The  number  of  elements  in  the  2"^  dimension  of  ARRAY. 
Error  flag  (=0  for  no  error,  =1  if  field  was  not  present). 


5.4.9  Subroutine  DAOUT 

This  routine  writes  the  field  as  specified  by  name  and  level.  If  IRET=0,  the  data  was  extracted 
successfully,  IRET=1,  the  field  or  level  was  invalid  or  not  present,  IRET=2,  there  was  an  error 
while  reading  the  file. 

Calling  Sequence:  subroutine  daout(lu,idir,rdir,cdir,name,rlev,array,ni,nx,ny,iret) 

Data  Declaration:  logical  lu 

real  rdir,  rlev,  array 
integer  idir,  ni,  nx,  ny,  iret 
char  name,  cdir 

Arguments:  See  DAINP 

(b)  Sequential  Files 

GCOM2D  uses  a  sequential  file  format  for  storage  of  output  data.  The  subroutines  associated  with 
sequential  file  creation  and  access  are  stored  in  STSUBS.FOR 

5.4.10  Subroutine  OPNINP 

This  routine  opens  a  sequential  file  and  reads  the  parameter  arrays. 

Calling  Sequence:  subroutine  opninp(lunit,fhame,title,iparam,rparam.nparam) 

Data  Declaration:  logical  iunit 

char  fnamc,  title 
integer  iparam,  nparam 
real  rparani 
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Arguments:  lunit 
fiiame 
title 
iparam 
rparam 
nparam 


Logical  unit  of  sequential  file. 

Character  variable  containing  file  name. 
80  character  title. 

Integer  file  header  information. 

Real  file. 

Dimension  of  IPARAM  and  RPARAM. 


5.4.11  Subroutine  SEQINP 


This  routine  reads  the  fields  fi'om  a  sequential  file. 

Calling  Sequence:  subroutine  seqmp(Iunit.varuii,nx,ny,k, title, hour, mmut.day, month, 


Data  Declaration:  logical  lunit 
real  var 

integer  ni,  nx,  ny,  k,  hour,  minut,  day,  month,  year 
char  title 


Arguments:  lunit 
var 
ni 
nx 
ny 
k 

title 

hour,  minut, 
day,  month,  year 


Logical  unit  number  of  the  sequential  file. 

Array  to  be  read  in. 

The  first  dimension  of  VAR. 

The  number  of  elements  in  the  1  dimension  of  VAR 
The  number  of  elements  in  the  2"^  dimension  of  VAR. 
The  vertical  level. 

80  character  title  for  labeling. 

The  time  to  which  the  input  array  corresponds. 


5.4.12  Subroutine  OPNOUT 


This  routine  opens  a  sequential  file  and  writes  the  parameter  arrays. 

Calling  Sequence:  subroutine  opnout(lunit,fhame,title,iparam,rparam,nparam) 

Data  Declaration:  logical  lunit 

real  rparam 

integer  iparam,  nparam 
char  fhame,  title 

Arguments:  Sec  OPNINP 
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5.4.13  Subroutine  SEQOUT 

This  routine  writes  the  fields  to  a  sequential  file. 

Calling  Sequence:  subroutine  seqout(lunit,var,ni,nx,ny,k, title, hour, minut, day, month, 

year) 

Data  Declaration:  logical  limit 

real  var,  title 

integer  ni,  nx,  ny,  k,  hour,  minut,  day,  month,  year 
Arguments:  See  SEQINP 
(c)  Time  Series  Files 

A  standard  ASCII  format  is  used  for  storing  station  information  in  both  the  atmospheric  and  ocean 
models.  Model  output  written  in  this  format  can  be  subsequently  plotted  using  the  display  program 
TSPLOT. 

5.4.14  Subroutine  TSDOUT 

This  subroutine  writes  a  standard  time  series  file.  The  five  standard  fields  stored  in  atmospheric  time 
series  files  are  station  pressure,  wind  speed,  wind  direction,  temperature  and  rainfall.  The  standard 
three  fields  stored  in  ocean  model  time  series  files  are  sea  surface  height,  current  speed,  current 
direction.  The  remaining  two  columns  of  data  are  nominally  saved  for  ocean  temperature  and 
salinity  \\  hich  is  not  relevant  in  the  current  application. 

Calling  Sequence:  subroutine  tsdout  (lu,stname,stlat,stlong,t2one,depth,ihour, 

iminut,iday,imonth,iyear,  spd,dir,z,temp,salin,nptsJnter,source) 

Data  Declaration:  real  stlat,  stlong,  tzone,  depth,  spd,  dir,  z,  temp,  salin 

integer  ihour,  iminut,  iyear,  iday,  imonth,  npts,  inter 
char  stname,  source 
logical  lu 
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Arguments; 


lu 

stname 

stlat 

stlong 

tzone 

depth 


Logical  unit  number  of  time  series  output  file. 

80  character  string  containing  station  name. 

Station  latitude. 

Station  longitude. 

Wind  time  zone,  in  real  hours  relative  to  GMT. 

Station  depth  above  (-ve)  or  below  (+ve)  mean  sea  level 
ihour,  iminut,  iday, 

imonth,  iyear  The  time  that  the  output  arrays  correspond  to. 
spd  Current  speed  array. 

Current  direction  array. 

Sea  level  array. 

Inactive. 

Inactive. 

Dimension  of  output  arrays. 

Timestep  of  output  in  minutes. 

Character  string  describing  source  of  time  series  output. 


spd 
dir 
z 

temp 

salin 

npts 

inter 

source 


5.4.15  Subroutine  TSDINP 

The  purpose  of  this  subroutine  is  to  read  a  standard  time  series  file. 

Calling  Sequence:  subroutine  tsdinp(lu,stname, stlat, stlong, tzone, depth, 
ihour,iminut,iday,imonth,iyear, 
spd, dir,z,temp,salin,nobs,npts,inter, source) 

Data  Declaration;  real  stlat,  stlong,  tzone,  depth,  spd,  dir,  z,  temp,  salin 

integer  ihour,  iminut,  iyear,  iday,  imonth,  nobs,  npts,  inter 
char  stname,  source 
logical  lu 

Arguments:  See  TSDOUT 


(d)  Topography  Files 

A  standard  ASCII  fomat  is  used  for  writing  the  topography  and  bathymetiy  of  a  model  urid 
Output  wntten  in  this  format  can  be  subsequently  plotted  using  tlie  display  program  TOPLOT. 

5.4.16  Subroutine  TOPOUT 

This  subroutine  \^tes  a  standard  topography  file.  The  units  are  in  meters  and  the  zero  elevation 
denotes  mean  sea  level.  Topography  is  positive  and  bathymetry  is  negative. 

Calling  Sequence:  subroutine  topout(lu,titleoproj,proj,height,ni,nx,ny) 
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Data  Declaration:  real  proj,  height 

integer  jproj,ni,  nx,  ny 
char  title 
logical  lu 

LxDgical  unit  number  of  input  file. 

80  character  string  for  description  of  file. 

Projection  parameter  flag. 

7  element  array  of  projection  parameters. 

Array  of  heights. 

The  first  dimension  of  HEIGHT. 

The  number  of  elements  in  the  1st  dimension  of  HEIGHT. 
The  number  of  elements  in  the  2"^  dimension  of  HEIGHT. 

5.4.17  Subroutine  FLDIN 


Arguments:  lu 
title 
jproj 
proj 
height 
ni 
nx 
ny 


The  purpose  of  this  subroutine  is  to  read  a  standard  topography  file. 

Calling  Sequence:  subroutine  fldin(lu,region  jproj, proj, array, ni,nx,ny) 

Data  Declaration:  real  proj,  array 

integer  jproj,  ni,  nx,  ny 
char  region 
logical  lu 


Arguments:  lu 

region 

jproj 

proj 

array 

ni 

nx 

ny 


Lx)gical  unit  number  of  input  file. 

80  character  string  for  description  of  file. 

Projection  parameter  flag . 

7  element  array  of  projection  parameters. 

Array  of  heights. 

The  first  dimension  of  ARRAY. 

The  number  of  elements  in  the  1st  dimension  of  ARRAY. 
The  number  of  elements  in  the  2"^  dimension  of  ARRAY. 


5.5  Projection  Routines 

In  the  PCTides  system,  boUi  MAPS  and  GCOM2/3D  make  use  of  a  standard  Cartesian  grid  system 
for  the  solution  of  the  equations  of  motion.  The  eurvature  of  the  earth  is  taken  into  aecount  by  the 
use  of  a  Lambert  eonformal  map  projection  that  consists  of  projecting  the  surface  of  the  earth  onto  a 
nght  circular  cone  and  flattening  it  onto  a  plane  surface.  This  proeess  results  in  a  map  scale  faetor 
m,  which  is  defined  by  the  ratio  of  tlie  distance  on  the  grid  to  the  corresponding  distanee  on  the 
earth’s  surface.  The  subroutines  associated  with  grid  projection  are  described  as  follows. 
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5.5.1  Subroutine  GRDPRJ 


This  routine  initializes  the  grid  projection.  A  standard  Lambert  conformal  map  projection  is  used 

projecuon, 

folbws-  ^  dimension  7  that  contains  the  map  projection  information  as 


Calling  Sequence:  subroutine  grdpij(jproj,proj) 

Data  Declaration:  real  proj 

integer  jproj 


Arguments:  jproj 
proj 

PROJ(l) 

PROJ(2) 

PROJ(3) 

PROJ(4) 

PROJ(5) 

PROJ(6) 

PROJ(7) 


Projection  parameter  flag . 

7  element  array  of  projection  parameters. 

STLATl :  First  standard  latitude. 

ANML:  Normal  longitude. 

STLAT2:  Second  standard  latitude. 

DU:  Number  of  grid-points  from  the  southwest  comer  of  erid  to 
ANML.  ^ 

TANL:  Latitude  that  forms  a  tangent  with  ANML  at  southernmost 
point  on  grid. 

NY:  Number  of  grid  rows. 

STDXY :  Standard  grid  spacing  in  kms. 


5.5.2  Subroutine  GRDLL 


The  purpose  of  this  subroutine  is  to  return  the  (i  j)  coordinates  for  a  specified  latitude  and  longitude. 

Calling  Sequence:  entry  grdll(jproj,rlat,rlong,ri,ij,theta) 

Data  Declaration:  real  rlat,  rlong,  ri,  ij,  theta 
integer  jproj 


Arguments:  jproj 
rlat 
rlong 
ri 

theta 


Projection  parameter  flag  (set  to  3  for  Lambert  conformal). 
Latitude  (supplied). 

Longitude  (supplied). 

The  i-coordinate  (returned). 

The  j -coordinate  (returned). 

Inactive. 
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5.5.3  Subroutine  GRDIJ 

The  purpose  of  this  subroutine  is  to  return  the  latitude  and  longitude  for  a  specified  (ij)  point. 

Calling  Sequence;  entry  grdij(jproj,rlat,rlong,ri,ij, theta) 

Data  Declaration:  real  rlat,  rlong,  ri,  tj,  theta 
integer  jproj 

Arguments:  See  GRDLL 


5.5.4  Subroutine  GRDLIM 

The  purpose  of  this  subroutine  is  to  calculate  the  latitude  and  longitude  limits  of  a  grid. 

Calling  Sequence:  subroutine  grdlim(jproj,nx,ny, slat, nlat,wlong,elong, region) 

Data  Declaration:  real  slat,  wlong,  elong 

integer  jproj,  nx,  ny,  nlat 
char  region 


Arguments;  jproj 

Projection  parameter  flag . 

nx 

Number  of  grid-points  in  x-direction. 

ny 

Number  of  grid-points  in  y-direction. 

slat 

Southern  latitude  limit  of  grid. 

nlat 

Northern  latitude  limit  of  grid. 

wlong 

Western  latitude  limit  of  grid. 

elong 

Eastern  latitude  limit  of  grid. 

region 

80  character  string  for  description  of  file. 

5.6  Date/Time  Routines 

Date  and  time  routines  are  used  to  convert  between  date  and  time  specified  on  input  files  and  an 
invariant  reference  time  used  by  the  model.  The  reference  time  used  by  the  PCTides  models  is  the 
P' January  1900. 

5. 6. 1  Subroutine  A  TIME 

The  purpose  of  this  subroutine  is  to  return  a  double  precision  number  of  hours  since  1*'  January 
1 900  for  a  supplied  hour  (hh),  minute  (mm),  day  (DD),  montli  (MM)  and  year  (YYYY). 

Calling  Sequence:  subroutine  atime(hour,minut,day,nionth,year,t) 
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Data  Declaration;  real  t 

integer  hour,  minut, 


Arguments:  hour 
minut 
day 
month 
year 
t 

5. 6. 2  Subroutine  ADA  TE 


day,  month,  year 

Hour  (supplied). 

Minute  (supplied). 

Day  (suppli^). 

Month  (supplied). 

Year  (supplied). 

Number  of  hours  since  origin. 


year  (YYYY)  when  supplied  with  a  double  precision  number  of  hours  since  January  1900. 

Calling  Sequence:  entry  adate(hour,minut,day,month,year,t) 

Data  Declaration:  real  t 

integer  hour,  minut,  day,  month,  year 

Arguments:  See  ATIME 

5, 6. 3  Subroutine  A  TIME2 


The  purpose  of  this  subroutine  is  to  return  a  double 


precision  number  of  hours  since  1"'  Januaiy 


1 900  when  supplied  with  the  date  in  (YYYYMMDD)  foimat  and  the  time  in  (hhmm)  format^ ' 
Calling  Sequence:  subroutine  atime2(date,time,t) 

Data  Declaration;  real  date,  time,  t 


Arguments:  date 
time 
t 

5.6.4  Subroutine  ADATE2 


YYYYMMDD. 

HHMM. 

Time  since  origin. 


The  puipose  of  this  subroutine  is  to  return  the  date  in  (YYYYMMHH)  format  and  lime  in  (hhmm) 
fonnat  \\  hen  supplied  with  a  double  precision  number  of  hours  since  1  January  1 900. 

Calling  Sequence:  entiy  adate2(date,timc,t) 

Data  Declaration:  real  date,  time,  t 

Arguments:  See  AT1ME2 


80 


PSI  Technical  Rq>ort  SSC-005-00 


PCTides  SDD 


5.7  Filtering  Routines 

It  is  common  practice  in  finite-difference  models  of  the  atmosphere  and  ocean  to  apply  numerieal 
filters  to  the  model  simulated  variables  (Haltiner  and  Williams,  1980).  This  is  to  reduce  numerieal 
errors  and  computational  instabilities  that  may  lead  to  a  spurious  growth  of  short  waves  that  may 
mask  the  actual  numerical  simulation  or,  in  severe  cases,  eould  eause  eatastrophie  instability.  In 
GCOM2D,  the  following  five-point  “low-pass”  filter  (Shapiro,  1970)  is  applied  to  the  model 
variables  in  each  grid  direction  approximately  onee  every  two  to  ten  physies  time-steps; 


where  0  represents  the  particular  model  variable  being  filtered  and  (p*,  the  filtered  variable.  The 
coefficients  a,  b,  and  c  are  -0.0625,  0.25  and  0.625,  respeetively.  The  filter  specifically  targets  two- 
grid-inten,  al  and  some  sub-grid-interval  waves,  damps  all  other  waves  except  for  the  infinitely  long 
wave.  Wa\'e  phases  remain  unchanged  by  the  filtering  process. 

5. 7. 1  Subroutine  HFILT 

HFILT  filters  the  entire  array  that  is  supplied  to  it. 

Calling  Sequence:  subroutine  hfilt(aiTay,ni,nx,ny,rl) 

Data  Declaration:  real  array,  rl 

integer  ni,  nx,  ny 

Array  of  values  to  be  filtered. 

The  first  dimension  of  ARRAY. 

The  number  elements  in  the  dimensions  of  ARRAY. 

The  number  elements  in  the  2"^  dimensions  of  ARRAY. 

Inactive. 

5. 7.2  Subroutine  HFIL T3 


Arguments:  array 
ni 
nx 
ny 
rl 


HF1LT3  is  used  for  filtering  ocean  model  fields  that  are  interrupted  by  irregular  coastlines.  Land 
points  are  defined  by  special  values. 

Calling  Sequence:  subroutine  hfilt3(array,ni,nx,ny,rl,spval) 

Data  Declaration:  real  array,  rl,  spval 
integer  ni,nx,  ny 
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Arguments:  array 
ni 
nx 
ny 
rl 

spval 


Array  of  values  to  be  filtered. 

The  first  dimension  of  ARRAY. 

The  number  elements  in  the  dimensions  of  ARRAY. 
The  number  elements  in  the  2^  dimensions  of  ARRAY. 
Inactive. 

special  value  used  to  mask  out  elements  in  an  array. 


5.8  Miscellaneous  MAPS  Routines 


There  are  four  subroutines  pertaining  to  the  solution  of  the  semi-implicit 
reside  in  the  AMSUBS  library. 


equations  in  MAPS  that 


5. 8. 1  Subroutine  MA  TINV 


This  routine  inverts  a  matrix  using  Gauss-Jordon  pivotal  elimination. 

Calling  Sequence:  subroutine  matinv(a,b,nmax,n,l,d,irror) 

Data  Declaration:  real  a,  b,  d 

integer  nmax,  n,  1,  irror 


NxN  matrix  to  be  inverted  on  entry,  inverse  on  exit. 

NxN  matrix  =  RHS  of  the  equation  on  entry,  solution  on  exit. 

The  first  dimension  of  the  arrays  A  and  B. 

The  number  of  elements  stored  in  A  and  B. 

If  ls>0,  solutions  are  given; 

L=0,  inverse  is  given; 

L<0,  both  are  given 
Determinant  of  matrix. 

=1  if  matrix  is  singular 
=0  if  matrix  is  inverted  successfully. 

5.8.2  Subroutine  OPTMUM 

The  purpose  of  this  subroutine  is  to  optimize  the  over-relaxation  coefficient  for  the  routine  LIEBH 
CalUng  Sequence:  subroutine  optmum(idimjdim,alfa,sigma) 

Data  Declaration:  real  alfa,  sigma 

integer  idim,  jdim 


Arguments:  a 
b 

nmax 

n 

1 


d 

irror 


Arguments;  idim 
jdim 
alfa 
sigma 


NX-1. 

NJ-1. 

Relaxation  coefficient  calculated  in  OPTMUM. 

Input  required  for  calculation  of  relaxation  coefficient. 
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5.8.3  Subroutine  EIGENP 


The  purpose  of  this  subroutine  is  to  find  the  eigenvalues  and  the  eigenvectors  of  a  real  general 
matrix  of  order  N. 

Calling  Sequence:  subroutine  eigenp(a, veer, veci,evr,evi,nm,n, indie) 

Data  Declaration:  real  a,  veer,  veci,  evr,  evi 
integer  nm,  n,  indie 

Matrix  of  N  rows  and  N  columns. 

NxN  matrix  of  real  eigenvectors  calculated  in  the  routine. 

NxN  matrix  of  imaginary  eigenvectors. 

Vector  of  N  dimension  containing  the  real  eigenvalues. 

Vector  of  N  dimension  containing  the  imaginary  eigenvalues. 

The  first  dimension  of  the  arrays  A,  VECR,  VECI. 

The  number  of  elements  stored  in  A,  VECR,  VECI,  EVR,  EVI. 

=0  when  no  eigenvalues  or  vectors  are  found 
=1  when  eigenvalues  are  found  but  eigenvectors  are  not  found 
=2  when  both  eigenvalues  and  eigenvectors  are  found. 

5.8.4  Subroutine  LIEBH 


Arguments:  a 

veer 

veci 

evr 

evi 

nm 

n 

indie 


The  purpose  of  this  subroutine  is  to  solve  the  Helmholtz  equations. 

Calling  Sequence:  subroutine  liebh(array,sol,emsqi,alf,itns,sigma,iI  jl,ax,ny,nz) 

Data  Declaration:  real  array,  sol,  emsqi,  alf,  sigma 
integer  ims,  il,  jl,  nx,  ny,  nz 


Arguments:  array 

sol 

emsqi 

alf 

ilns 

sigma 

il 

Jl 

nx 

ny 

nz 


ILxJLxNZ  array  of  input  to  subroutine. 

ILxJLxNZ  array  containing  solution  of  Helmholtz  equations. 
Inverse  of  map  factor  squared. 

Vector  of  NZ  relaxation  coefficients. 

Vector  of  No.  of  iterations  required  for  solution  at  each  level. 
Vector  of  NZ  eigenvalues  multiplied  by  As'. 

The  first  dimension  of  ARRAY  and  SOL. 

The  second  dimension  of  ARRAY  and  SOL. 

No.  of  elements  in  first  dimension  of  ARRAY  and  SOL. 

No.  ot  elements  in  second  dimension  of  .A.RRAY  and  SOL. 
The  third  dimension  of  ARRAY  and  SOL.  No.  of  vert,  levels. 
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5.9  Miscellaneous  GCOM  Routines 

There  are  several  subroutines  called  by  GCOM2/3D  that  are  stored  in  the  library  STSUBS.FOR 

5.9.1  Subroutine  WEIGHT 


The  purpose  of  this  subroutine  is  to  calculate  the  boundary  weight  factors  for  nesting  and 
application  of  sponge  layers.  ^ 

Calling  Sequence:  subroutine  weight(npts, linear, wgt,ni,nx,ny) 

Data  Declaration:  real  wgt 

integer  npts,  linear,  ni,  nx,  ny 

Number  of  nesting  points. 

Interpolation  flag  (internally  set  to  1). 

Array  of  nesting  weights. 

The  first  dimension  of  depz. 

The  number  of  elements  stored  in  the  dimension  of  depz. 

The  number  of  elements  stored  in  the  dimension  of  depz. 

5. 9. 2  Subroutine  SETWGT 


Arguments:  npts 
linear 
wgt 
ni 
nx 
ny 


The  purpose  of  this  subroutine  is  to  set  the  nesting  weights. 

Calling  Sequence:  subroutine  setwgt(wgt,depz,ni,nx,ny,npts) 

Data  Declaration:  real  wgt,  depz 

integer  ni,  nx,  ny,  npts 


Arguments:  wgt 

depz 

ni 

nx 

ny 

npts 


Array  of  nesting  weights. 

Array  of  water  depths  relative  to  mean  sea  level. 
The  first  dimension  of  depz  and  wgt. 

The  number  of  grid  points  in  x-direction. 

The  number  of  grid  points  in  y-direction. 
Number  of  nesting  points. 


system  clock. 


5. 9.3  Subroutine  START 

The  purpose  ol  this  subroutine  is  to  calculate  the  start  time  Ifom  the 
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5.9.4  Subroutine  THPRED 

The  purpose  of  this  subroutine  is  to  read  astronomieal  constants  for  each  tidal  constituent  and 
calculate  resulting  tidal  height. 

Calling  Sequence:  subroutine  thpred(lstc,con,ntcmax,ntcs,anip,phase,chp,ch,rlat, 
tdhgt,ni,nx,ny, hour, minut, day, month, year, dt,tzone) 


Data  Declaration:  logical  Istc 

real  con,  amp,  phase,  chp,  eh,  rlat,  tdhgt,  dt,  tzone 

integer  ntemax,  ntes,  ni,  nx,  ny,  hour,  minut,  day,  month,  year 


Arguments:  Istc 

Logical  unit  which  reads  the  standard  tidal  constituent  data  from 
STCONTHP.DAT. 

con 

Tidal  constituents  array,  CON(NTCMAX,NOBS). 

ntemax 

Maximum  number  of  tidal  constituents. 

tes 

Number  of  tidal  constituents  available  at  each  tidal  station. 

amp 

Tidal  amplitudes  read  in  from  M2.DAT,  S2.DAT  etc. 

phase 

Tidal  phases  read  in  from  M2.DAT  etc. 

chp 

Working  array  dimensioned  (NTCMAX,NX,NY). 

ch 

As  above. 

rlat 

Station  latitude. 

tdhgt 

Predicted  tidal  height. 

ni 

First  dimension  of  AMP,  PHASE  and  TDHGT. 

nx 

The  number  of  elements  in  the  1^  dimension  of  AMP,  etc. 

ny 

hour,  minut, 
day,  month. 

The  number  of  elements  in  the  2™*  dimension  of  AMP,  etc. 

year 

Time  of  required  tidal  height  prediction. 

dt 

Model  time-step. 

tzone 

Wind  time  zone,  in  real  hours  relative  to  GMT. 

5.10  Data  Flics  for  GCOM2D/3D  and  MAPS 


In  order  tor  MAPS  and  GCOM2/3D  to  be  globally  relocatable  there  needs  to  be  global 
climatological  and  environmental  data  files  to  establish  basic  grid  and  model  fields  and 
parameters.  These  data  files  are  summarized  below. 


5. 10. 1  Topographic/Bathymetric  Data 

Topographic/bathymctric  data  arc  stored  on  latitudc/longitude  grids  in  a  series  of  direct  access 
files,  all  with  the  file  extension  “.DA”.  All  values  arc  .stored  as  1NTEGER*2  and  the  direct 
access  file  structure  is  as  follows: 
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Record  1 ;  Southern  latitude  limit 
Record  2;  Northern  latitude  limit 
Record  3;  Western  longitude  limit  (0-360) 
Record  4;  Eastern  longitude  limit.  (0-360) 
Record  5:  Resolution  in  minutes  x  1000 
Record  6;  Topography  -  column  1,  row  1 
Record  7:  Topography  -  column  2,  row  1 
Etc. 

5.10.1.1  Subroutine  GETOP 


This  routine  returns  the  topography  VALUE  at  the  point  (RLAT,  RLONG)  from  the  hiehe<;t 
resolution  .DA”  file  eovering  the  point.  GETOP  is  in  the  library  file  DASUBS.FOR. 

Calling  Sequence;  getop(rlat,rlong, value) 

Data  Declaration:  real  rlat,  rlong,  value 


Arguments: 


rlat 

rlong 

value 


Station  latitude. 
Station  Longitude. 
Topography  value. 


5.10.2  Climatological  Data 


MAPS  also  requires  some  global  climatological  data  for  roughness  length,  albedo  and  soil 
motsture.  Roughness  length,  albedo  and  monthly  soil  moisture  values  are  stored  in  INTEGER'*2 

?r'^OIl‘\l‘vPRn“‘  GLROUGHLDA,  GLALBEDO.DA  and  GLSOILMV  JAN 

GLSOILMV.FEB  etc.  The  direct  access  file  structure  is  as  follows: 


Record  1 :  Southern  latitude  limit 
Record  2;  Northern  latitude  limit 
Record  3:  Western  longitude  limit  (0-360) 
Record  4;  Eastern  longitude  limit.  (0-360) 
Record  5:  Resolution  in  minutes 
Record  6:  Factor  to  multiply  data  by 
Record  7;  Value -column  1.  row  1 
Record  8:  Value  —  column  2,  row  I 
Etc. 


5.10.2.1  S  uhroutine  DA  2D  RE  A  D 


mLATRLONrf  nA%nnpTr''-  T  ''''‘•UE  at  the  point 

(KLA 1 , RLONG).  DA2DREAD  is  in  the  library  file  DASUBS.FOR. 
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Calling  Sequence:  da2dread  (lu,fhame,rlat,rlong, value) 

Data  Declaration;  logical  lu 
char  fhame 
real  rlat,  rlong,  value 

Arguments:  lu  Lx)gical  unit  number  of  the  model  run  log  file, 

fhame  Character  variable  containing  file  name, 

rlat  Station  latitude, 

rlong  Station  Longitude, 

value  Value. 

The  program  DAF1LE.EXE  provides  the  ability  to  convert  these  files  between  an  ASCII  format 
and  there  direct  access  binary  format.  The  program  recognizes  the  extension  “.DA”  as  the  direct 
access  file  and  the  extension  “.ASC”  as  the  ASCII  file. 

5.10.3  Tidal  Data 

GCOM2/3D  also  require  global  tidal  files  from  the  FES95. 1/2.1  model  to  derive  boundary  tidal 
forcing.  These  30  minute  resolution  data  are  stored  in  the  ASCII  files  GLBM2.30M, 
GLBS2.30M,  etc.  These  routines  are  in  the  library  file  STSUBS.FOR. 

5.10.3.1  Subroutine  GETIDE 

This  routine  reads  these  ASCII  files  and  returns  the  amplitude,  AMP,  phase,  PHASE,  and  period, 
PERIOD,  of  the  tidal  constituent,  CON,  at  the  point  (RLAT,RLONG),  in  time  zone,  TZONE. 
GETIDE  is  in  the  library  file  DASUBS.FOR. 

Calling  Sequence:  getide(rlat, rlong, con, amp, phase,period,tzone) 

Data  Declaration:  real  rlat,  rlong,  con,  amp,  phase,  period,  tzone 

Arguments:  rlat  Station  latitude. 

rlong  Station  Longitude, 

con  Tidal  constituent, 

amp  Amplitude, 

phase  Phase, 

period  Period, 

tzone  Time  zone. 


5.10.3.2  Subroutine  TIDEOUT 

The  purpose  of  this  routine  is  to  write  an  ASCII  tidal  file  on  logical  unit,  LU,  for  the  region. 
REGION,  projeetion  parameters,  JPROJ,  PROJ(7),  tidal  eonstituent,  CON,  with  period. 
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PERIOD,  in  time  zone,  TZONE,  and  amplitude,  AMP,  and  phase,  PHASE 
dimensions,  NX,  NY.  ’ 


arrays  with 


CaUing  Sequence:  tidout(lu,region,jproj,proj,con,period,tzone,amp,phase,ni,nx,ny) 

Data  Declaration:  logical  lu 

char  region 

integer  jproj,  ni,  nx,  ny 

real  proj,  con,  period,  tzone,  amp,  phase 


Arguments: 


lu 

Logical  unit  number  of  the  model  run  log  file. 

region 

80  character  title  fi-om  the  TOPOG.DAT  file. 

jproj 

Flag  defining  Lambert-conformal  projection. 

proj 

Array  of  Lambert-conformal  map  projections. 

con 

Tidal  constituent. 

period 

Period. 

tzone 

Time  zone. 

amp 

Amplitude. 

phase 

Phase. 

ni 

Maximum  horizontal  dimension  (x -direction). 

nx 

The  number  of  grid  points  in  x-direction. 

ny 

The  number  of  grid  points  in  y-direction. 

5.10.3.3  Subroutine  TIDEIN 


The  purpose  of  this  subroutine  is  to  read  an  ASCII  tidal  file. 

Calling  Sequence:  tidein(lu,region,jproj,proj,con,period,tzone,amp,phase,ni,nx,ny) 

Data  Declaration:  logical  lu 

real  region,proj,  con,  period,  tzone,  amp,  phase,  ni,  nx,ny 
integer  jproj 


Arguments: 


lu 

Logical  unit  number  of  the  model  run  log  file. 

region 

80  character  title  Ifom  the  TOPOG.DAT  file. 

jproj 

Flag  defining  Lambert-conformal  projection. 

proj 

Array  of  Lambcrt-confonnal  map  projections. 

con 

Tidal  constituent. 

period 

Period. 

tzone 

Time  zone. 

amp 

Amplitude. 

phase 

Phase. 

ni 

Maximum  horizontal  dimension  (x-direction). 

nx 

The  number  of  grid  points  in  x-direction. 

ny 

The  number  of  grid  points  in  y-direction. 

88 


PSI  Technical  Report  SSC-005-00 


PCTides  SDD 


6.0  REQUIREMENTS  TRACEABILITY 

The  Software  Test  Description  (STD)  includes  eight  test  runs  that  encompass  all  CSCI  routines 
(SDD)  and  requirements  (SRS). 


SRS  paragraph  numbers: 


SDD 

Paragraph 

numbers: 

3.1.1 

3.1.2 

3.1.3 

3.1.4 

3.1.5 

vrre 

3.1.2A 

3.1.3A 

3.1.4A 

5.1 

X 

X 

X 

5.1.2 

X 

5.1.3 

X 

5.1.4 

X 

X 

X 

X 

5.2 

X 

5.2.2 

X 

5.2.3 

X 

X 

5.2.4 

X 

X 

X 

X 

X 

5.3 

X 

5.3.2 

X 

5.3.3 

X 

5.4 

X 

X 

5.5 

X 

5.6 

X 

5.7 

X 

5.8 

X 

5.9 

X 

5.10.1 

X 

X 

5.10.2 

X 

X 

5.10.3 

X 

89 


PSI  Technical  Report  SSC-005-00 


PCTides  SDD 


7.0  NOTES 


7.1  Acronyms  and  Abbreviations 


ASA 

CSCI 

CSC 

FES 

GCOM2D 

GCOM3D 

GMT 

MAPS 

NOGAPS 

NRL 

OAML 

PC 

PCTides 

PSI 

SDD 

SRS 

SSC 

STD 

UNIX 


Applied  Sciences  Associates 

Computer  Software  Configuration  Item 

Computer  Software  Component 

Finite  Element  Solution 

Coastal  Ocean  Model  2-D 

Coastal  Oeean  Model  3-D 

Generic  Mapping  Tool 

Mesoscale  Atmospheric  Prediction  System 

Navy  Operational  global  Atmospheric  Prediction  System 

Naval  Researeh  Laboratory 

Oceanographic  and  Atmospheric  Master  Library 

Personal  Computer 

Globally  Relocatable  Navy  Tide/Atmospheric  Modeling  System 

Planning  Systems,  Incorporated 

Software  Design  Description 

Software  Requirements  Specification 

Stennis  Space  Center 

Software  Test  Description 

Workstation  Operating  System 
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8.0  Appendix  I.  FORTRAN  Common  Blocks 

APPENDIX  I.  FORTRAN  COMMON  BLOCKS .  91 

CSC  Maps  COMMON  BLOCKS .  92 

COMMON/  DOUBLE/ .  92 

COMMON/ THREED/ .  92 

COMMON/ TWOD/ .  93 

COMMON/ ONEDR/ .  94 

COMMON/  INTEGR/ .  95 

COMMON/  ONEDC/ .  95 

CSC  GCOM2D  COMMON  BLOCKS .  96 

COMMON/ ARRAYS/ .  96 

CSC  GCOM3D  COMMON  BLOCKS .  97 

COMMON/  THREED/ .  97 

COMMON/ TWOD/ .  97 
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CSC  Maps  COMMON  Blocks 


COMMON  /  DOUBLE/ 

precision  arrays  (both  two-dimensional  and  three- 


Variable 

-Type 

Description 

sol  (m,  nj,  nlc),  array  (ni,  nj,  nk) 

Real 

Arrays  used  for  solution  of  Helmholtz  equations;  array 
contains  the  forcing  function. 

emsqi  (ni,nj) 

Real 

Inverse  of  the  square  of  the  map  factor 

psm  (m,  nj),  ps(m,nj), 
_psp(ni,nj) 

Real 

Surface  pressure  at  the  pre\aous,  current,  and  next  time 
steps. 

amv(nk,  nk),  fcainv  (nk,  nk), 
feainv  (nk,nk),  h  (nk,nk),  finat 
(nk,nk).  veer  (nk,nk),  vecrl 
(nk,  nk),  hinv  (nk,nk),  tinv 
(nk,nk),  tmat  (nk,nk),  tlag 
(nk,nk).  ulag  (nk,nk) _ 

Real 

Double  precision  matrices  used  in  the  semi-implicit 
solution  procedure. 

vecqq  (nk),  bam  (nk),  qvec 
(nk),  alfcor  (nk),  dq  (nk),  alf 
(nk),  e\Teal  (nk) 

Real 

Double  precision  vectors  used  in  the  semi-implicit 
solution  procedure. 

sigma  (nk) 

Real 

Array  for  calculating  relaxation  coefficientc 

dqmqa>  (nk) 

Real 

Convective  adjustment  array. 

wts  (nkp  1 ) 

Real 

Array  used  in  semi-implicit  solution  procedure 

ds,  dsi,  dsi2,  dssq 

Real 

The  horizontal  grid  spacing;  A  s,  and  related  parameters' 
As\{2  Asy^A  s^. 

bet65,dmiax,  alftim 

Real 

Constants  for  semi-implicit  scheme 

COMMON  /  THREED/ 

This  common  block  contains  all  three-dimensional  arrays. 


Variable 

-Type 

Description 

rmm  (nk.ni,nj),  nn  (nk,ni,nj), 
imp  (nk.ni,nj) 

Real 

Arrays  for  mixing  ratio  at  pre\aous,  current  and  next  time 
step. 

urn  (nk.ni.nj),  u  (nk,ni,nj), 
up(nk,ni.nj) 

Real 

Arrays  for  U-component  of  \  elocity  at  previous,  current, 
and  next  time  step. 

vm  (nk.ni.nj),  v  (nk,ni,nj),  vp 
(nk,ni,nj) 

Real 

Arrays  for  U-component  of  \  elocity  at  previous,  current,  : 
and  next  time  step.  ’  1 

phi(nk,ni.nj),  oincga(nk,ni,ni) 

Real 

Arrays  for  geopotential  height  and  vertical  velocity  i 

tm(nk,ni.nj),  t(nk,ni,nj), 
tp(nk,ni.nj) 

Real 

Arrays  for  temperature  at  pre\ious,  current  and  next  time  i 
step. 

sigdot(nk.ni,nj),  (dot(nk,ni,nj), 
qdot(nk.ni.nj) 

Real 

Arrays  for  cr ,  7\  and  /? , 
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COMMON  /  TWOD/ 

This  common  block  contains  all  two-dimensional  arrays. 


Variable 

Type 

Description 

trainm(ni.nj),  train(ni,nj), 
trainp(ni.nj) 

Real 

Arrays  for  total  rainfall  at  the  previous,  current,  and  next 
time  step  (convective  +  large  scale). 

tsm(ni,nj).  ts(ni,nj),  tsp(ni,nj) 

Real 

Arrays  for  soil  temperature  at  the  first  soil  layer  at  the 
previous,  current  and  next  time  step. 

tdm(ni,nj),  td(ni,nj),  tdp(ni,nj) 

Real 

Arrays  for  temperature  of  the  second  soil  layer  at  the 
previous,  current,  and  next  time  step. 

wsm(ni,nj),  ws(ni,nj), 
wsp(ni,nj) 

Real 

Arrays  for  soil  moisture  at  the  first  soil  layer  at  the 
previous,  current,  and  next  time  step. 

wdm(ni,nj),  wd(ni,nj), 
wdp(ni,nj) 

Real 

Arrays  for  soil  moisture  at  the  second  soil  layer  at  the 
previous,  current,  and  next  time  step. 

zs(ni,nj),  corp(ni,nj) 

Real 

Arrays  of  terrain  height  and  Coriolis  parameter. 

em(ni,nj),  emu(ni,nj), 
emv(ni,nj) 

Real 

Arrays  of  map  factor  at  the  temperature,  and  U  and  V 
(staggered)  grid  points. 

albd(ni,nj).  ruf(ni,nj), 
stmp(ni,nj).  smos(ni,nj) 

Real 

Arrays  for  surface  albedo,  roughness  length,  sub-surface 
temperature  and  sub-surface  moisture.  All  are  specified 
in  the  pre-processor  and  input  to  CSC  Maps. 

dummyl(ni.nj), 
dummy2(ni.nj), 
dummy3(ni,nj),  dummy4(ni,ni) 

Real 

Dummy  arrays  used  for  extra  calculations  in  various 
subroutines. 

qlat(ni,nj),  qlong(ni,nj) 

Real 

Arrays  for  latitude  and  longimde  of  each  grid  point. 

crain(ni,nj).  lrain(ni,nj) 

Real 

Arrays  for  convective  and  large  scale  rainfall. 

wgt(ni,nj).  w2(ni,nj) 

Real 

Weights  used  for  nesting  scheme.  w2=(1.0-wgt) 

hlatn(ni,ni).  hsens(ni,nj) 

Real 

Arrays  for  latent  and  sensible  heat. 
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COMMON  /  ONEDR/ 

This  common  block  contains  real  arrays  and  variables. 


Variable 

-Type 

Description 

tbar  (nk) 

Real 

Average  temperature  for  each  sigma  level.  Set  in 

NESTIN. 

q(nk) 

Real 

Sigma  levels.  Set  in  CONST. 

rtbar  (nk) 

Real 

Real  variable  for  tbar. 

teool  (nk) 

Real 

Newtonian  eooling  to  space  (set  in  PARAIU) 

qph  (nk),  qipk(nkDl) 

Real 

Half  sigma  levels,  storage  arrav 

bet  (nk),  dtodq  (nk),  gama  (nk) 

Real 

stlat  (nout),  stlong  (nout) 
stni  (nout),  stnj  (nout) 

Real 

latitudes  and  longitudes  and  corresponding  real  grid 
point  values  of  station  variables. 

istn  (nout).  jstn  (nout) 

Real 

Integer  grid  coordinates  of  output  stations 

philog  (nout) 

Real 

Grid  orientation  at  output  stations. 

proj  (7) 

Real 

Array  of  projection  parameters  defining  grid 

rdir  (nij/3) 

Real 

Array  of  real  numbers  stored  n  the  directoiy  header  of 
direct  access  input  and  output  files 

tfilt,  tfiltm 

Real 

Filter  coefficients  used  for  Asselin  Filter 

soldec 

Real 

Solar  declination  angle. 

hrtime 

Real 

Real  time  in  hours. 

rhf 

Real 

Relative  humidity  at  which  large  scale  rain  trieeers 

diffvM 

Real 

Honzontal  diffusion  coefficient  (not  nsed) 

hours,  drfacl,  dtop,  ffac, 
rmdfwt 

Real 

Run  time  in  hours,  time  step  reduction  factor,  default  top 
of  atmosphere,  filter  factor,  moisture  diffusion 
coefficient,  respectively. 

precp,  precta,  cks,  eke,  pe, 
psbar,  trhat.  vromg 

Real 

Precipitation  variables,  map  factor  storage  variable, 
kinetic  and  potential  energy  variables,  respectively 

cp,  g,  r  roncp 

Real 

Specific  heat  of  dry  air,  gravitational  constant,  gas 
constant,  R/Cp. 

dl,  d2.  d3 

Real 

Dcjpth  of  soil  layers  used  in  surface  temperature  and 
moisture  prediction  scheme. 

wwl,  \\^v2 

Real 

Soil  temperature  calculation  constants  for  second  layer 
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COMMON  /  INTEGR/ 

This  is  a  common  block  for  integer  arrays  and  variables. 


Variable 

Type 

Description 

land  (ni.nj) 

Integer 

Land/sea  mask  (1  for  land,  0  for  sea). 

itns  (nk) 

Integer 

Dummy  array  for  Helmholtz  solver. 

nstep 

Integer 

Integer  timestep,  initially  zero  in  P  ARAMS  and 
incremented  by  1  at  end  of  time-stepping  loop. 

nstepl,  nsiep2 

Integer 

First  and  last  time  steps. 

irst 

Integer 

Restart  time  step. 

nestep 

Integer 

Next  time  step  at  which  nesting  values  are  read. 

nphys 

Integer 

Flag  determining  how  frequently  the  physics  routine  is 
called. 

nest 

Integer 

Flag  for  nesting  scheme  (0=  open  boundaries  with 
radiation  condition,  1=  nesting  in  tendencies,  2=  nesting 
in  absolute  values. 

nsthrs,  nhours,  inthrs,  nb, 
modes 

Integer 

Nesting  time  in  hours,  total  run  time  in  hours,  nesting 
interval  in  hours,  number  of  nesting  grid  points,  resp. 

idate,  itime.  kdatc,  ktime, 
nstime,  ilog,  logi,  logj,  iplot 

Integer 

Analysis  date,  analysis  time,  forecast  date,  forecast  time, 
nesting  time,  resp. 

nx,  ny,  nz 

Integer 

Actual  model  grid  dimensions  in  the  x,  y,  and  z 
direction. 

nxy,  nzm  1 ,  nzm2,  nzm3, 
nzm4,  nzpl 

Integer 

Total  number  of  horizontal  grid  points,  number  of  sigma 
levels  minus  1 ,2,3,4  and  plus  1 . 

Iphys,  Isurf,  Iconls,  Iconcu, 
Iverdf,  Ihdiflf,  Ishcon 

Integer 

Flags  for  running  physics  scheme,  the  surface  scheme, 
large  scale  rain,  Kuo  scheme,  vertical  diffusion, 
horizontal  diffusion  and  shallow  convection,  resp. 

nstns 

Integer 

Number  of  stations  for  time  series  output. 

kdtop 

Integer 

Top  level  for  application  of  vertical  diffusion. 

ilaps 

Integer 

Integer  lapse  time. 

jproj 

Integer 

Flag  representing  the  required  map  projection 
(3=Lambert  confonnal,  9=  mercator). 

idir  (nij/3 ) 

Integer 

Array  storing  all  integer  information  for  the  directory  of 
a  direct  access  file. 

COMMON /ONEDC/ 

This  is  a  common  block  for  character  arrays. 


Variable 

Type 

Description 

edir  (nij/3) 

Char 

Array  storing  all  character  information  for  tlie  directory 
of  a  direct  access  file. 

stnamc  (nout) 

Char 

Array  storing  place  names  for  time  series  output. 
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CSC  GCOM2D  COMMON  Blocks 


COMMON/ ARRAYS/ 

This  is  the  main  common  block  for  GCOM2D.  It  canies  all  two-dimensional  arrays. 


Variable 

Type 

Description 

twx  (mmj),  twy(ni,ni).  tDafni.nj') 

Real 

Wind  stress  and  surtace  Drcssure  fendenriVc 

tz(ni,nj) 

Real 

Sea  level  nesting  tendencies. 

up  (0;ni,  0;nj),  vp  (0:ni,  0:nj),  zp 
(ni,nj) 

Real 

Velocity  components  and  sea  level  perturbation  at 
new  time  step. 

uo(0:ni,  0:nj),  vo(0:ni,  0:nj), 
zo(ni,nj) 

Real 

Velocity  components  and  sea  level  perturbation  at  old~ 
time  step. 

ubar(0;m,  0;nj),  vbarrO:ni.  O.-nj') 

Real 

.  U  at  V  grid  points  and  V  at  U  grid  ooints  ~ 

zres(m,nj),  zcg(ni,nj) 

Real 

Nesting  arrays. 

wspeed(ni,nj),  wdir(ni,nj), 
j)a(ni,nj) 

Real 

Wind  speed,  wind  direction  and  atmospheric  pressure. 

wx(0;m,  0:nj),  wy(0:ni,  0:ni) 

Real 

Wind  stress  components. 

pfz(m,nj),  pfu(0;ni,  0;nj),  pfV 
(0:ni,  0:ni) 

Real 

Map  factors  at  U,  V  and  C  grid  points. 

depz(mmj),  depu(0;ni,  0:nj), 
depv(0;ni,  Omj) 

Real 

Depth  at  U,  V,  and  ^  grid  points. 

height(mmj),  spfeat(ni,nj), 
topog2(nimj) 

Real 

Topography  arrays. 

cf(ni,nj).  l\'(ni,ni) 

Real 

Coefficient  of  bottom  fiiction,  Coriolis  force 

tdhgt(nimj),  e(ni,nj),  wgt(ni,ni) 

Real 

Tidal  height,  total  energy,  and  nestine  weipbts  rpcp 

dummyO(ni,nj),  dummyl(ni,nj), 
dummy2(ni,nj),  dummy3(0;ni, 
0;nj),  dummy4(0;ni,  0:nj) 

Real 

Working  arrays. 

gravl(ni.nj) 

Real 

Storage  array. 

iminz  (nj).  imaAz(nj), 
iminu(0;ni.  O.nj),  imaxu(0;ni, 
0;nj),  iminv(0:ni,  0:nj), 
imaxv(0;ni.  0;nj) 

Integer 

The  first  sea  point  on  west  of  gnd  and  last  sea  point 
on  east  of  grid  on  each  row  at  z,  u,  and  v  points. 

jproj 

Integer 

Projection  parameter  flag. 

proj  (7) 

Real 

Projection  parameter  arraw 

Hog 

Logical 

Logical  unit  of  log  file 

depmin,  sidx>m,  g,  hOcbf 

Real 

Minimum  bathymetrie  depth,  standard  grid  spacing 
(m),  acceleration  due  to  gravity,  default  \  alue  of 
bottom  friction,  respecti\  elv. 
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CSC  GCOM3D  COMMON  Blocks 


COMMON/ THREED/ 

This  is  the  main  common  block  for  GCOM3D.  It  carries  all  two-dimensional  arrays. 


Variable 

Type 

Description 

wp(nk,ni,nj),  wo(nk,ni,nj), 
wm(nk,ni,nj) 

Real 

Vertical  velocity  at  next,  present,  and  last  time  step. 

rhop  (nk,ni,nj) 

Real 

Density. 

cevu(nk,nimj),  cew(nk,ni,nj) 

Real 

Coefficient  of  eddy  viscosity  at  U  and  V  points. 

COMMON/ TWOD/ 

This  is  the  main  common  block  for  GCOMD.  It  carries  all  two-dimensional  arrays. 


Variable 

Type 

Description 

kmaxz  (0:ni,  0:nj),  kmaxu(0:ni, 
0:nj),  kmaxv  0:ni,  0:nj) 

Real 

Number  of  vertical  levels  at  each  Z,  U,  and  V  grid  point. 

bdepu(ni,nj),  bdepv(ni,nj), 
bdepz(ni,ni) 

Real 

Thickness  of  bottom  layer  at  U,  V,  and  Z  grid  points. 

h(nk) 

Real 

Thickness  of  each  vertical  la^’er. 

depk  (0;nk) 

Real 

Depth  of  each  vertical  layer. 

dzbz(ni,nj).  dzbu(ni,nj), 
dzbv(ni,nj) 

Real 

Distance  between  mid  point  of  bottom  layer  and  mid 
point  of  second  bottom  layer  at  Z,  U,  and  V  points. 
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